17#include "MantidHistogramData/HistogramE.h"
18#include "MantidHistogramData/HistogramY.h"
38using namespace Kernel;
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() {}
49const std::string NormaliseByPeakArea::name()
const {
return "NormaliseByPeakArea"; }
52int NormaliseByPeakArea::version()
const {
return 1; }
55const std::string NormaliseByPeakArea::category()
const {
return "CorrectionFunctions\\NormalisationCorrections"; }
62void NormaliseByPeakArea::init() {
63 auto wsValidator = std::make_shared<CompositeValidator>();
68 "An input workspace.");
70 auto mustBePositive = std::make_shared<BoundedValidator<double>>();
71 mustBePositive->setLower(0.0);
72 mustBePositive->setLowerExclusive(
true);
73 declareProperty(
"Mass", -1.0, mustBePositive,
"The mass, in AMU, defining the recoil peak to fit");
75 "If true all spectra on the Y-space, fitted & symmetrised workspaces "
76 "are summed in quadrature to produce the final result");
79 "Input workspace normalised by the fitted peak area");
81 "Input workspace converted to units of Y-space");
83 "Output from fit of the single mass peakin y-space. The output units are "
84 "in momentum (A^-1)");
86 "The input data symmetrised about Y=0. The output units are in momentum "
93void NormaliseByPeakArea::exec() {
98 const auto nhist =
static_cast<int64_t
>(yspaceIn->getNumberHistograms());
99 const auto nreports =
static_cast<int64_t
>(yspaceIn->getNumberHistograms() +
101 m_progress = std::make_unique<API::Progress>(
this, 0.10, 1.0, nreports);
103 for (int64_t i = 0; i < nhist; ++i) {
107 m_yspaceWS->setSharedX(i, yspaceIn->sharedX(i));
108 m_fittedWS->setSharedX(i, yspaceIn->sharedX(i));
112 double peakArea =
fitToMassPeak(yspaceIn,
static_cast<size_t>(i));
130void NormaliseByPeakArea::retrieveInputs() {
143 const size_t nhist =
m_sumResults ? 1 : yspaceIn->getNumberHistograms();
151 constexpr double high(1e6);
153 m_yspaceWS->setSharedX(0, yspaceIn->sharedX(0));
154 m_fittedWS->setSharedX(0, yspaceIn->sharedX(0));
171 auto xLabel = std::make_shared<Units::Label>(
"Momentum",
"A^-1");
186 alg->setProperty(
"InputWorkspace",
m_inputWS);
187 alg->setProperty(
"Mass",
m_mass);
194 double xmin(0.0), xmax(0.0);
195 tofInY->getXMinMax(xmin, xmax);
196 std::vector<double> params(3);
202 alg->setProperty(
"InputWorkspace", tofInY);
203 alg->setProperty(
"Params", params);
205 return alg->getProperty(
"OutputWorkspace");
219 func->setAttributeValue(
"Mass",
m_mass);
220 func->setAttributeValue(
"WorkspaceIndex",
static_cast<int>(
index));
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]);
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"
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);
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() <<
": "
263void NormaliseByPeakArea::normaliseTOFData(
const double area,
const size_t index) {
277 assert(yValues->rawData().size() == eValues->rawData().size());
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();
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)
293 double err = 1.0 / (accumEj * accumEj) + 1.0 / (rhsEj * rhsEj);
294 accumY[j] = accumYj / (accumEj * accumEj) + rhsYj / (rhsEj * rhsEj);
296 accumE[j] = 1.0 / sqrt(err);
299 accumWS->setSharedY(
index, yValues);
300 accumWS->setSharedE(
index, eValues);
307void NormaliseByPeakArea::symmetriseYSpace() {
317 const double dy = 0.1;
319 const auto nhist =
static_cast<int64_t
>(
m_symmetrisedWS->getNumberHistograms());
322 for (int64_t i = 0; i < nhist; ++i) {
329 for (
size_t j = 0; j < npts; ++j) {
330 const double ein = eyspace[j];
331 const double absXj =
fabs(xsym[j]);
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) {
341 double invE2 = 1 / (ek * ek);
344 double wt = (1 / (eout * eout)) + invE2;
#define DECLARE_ALGORITHM(classname)
IPeaksWorkspace_sptr workspace
std::map< DeltaEMode::Type, std::string > index
#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.
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
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.
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 m_inputWS
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.
double m_mass
The input mass in AMU.
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.
bool is(int level) const
Returns true if at least the given log level is set.
void information(const std::string &msg)
Logs at information level.
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.
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.
@ Input
An input workspace.
@ Output
An output workspace.