26const
std::
string FindEPP::name()
const {
return "FindEPP"; }
36 return "Performs Gaussian fits over each spectrum to find the Elastic Peak "
45 "An input workspace.");
47 "An output workspace.");
58 auto numberspectra =
static_cast<int64_t
>(
m_inWS->getNumberHistograms());
76 auto spectrum =
static_cast<size_t>(
index);
77 m_outWS->cell<
int>(spectrum, 0) =
static_cast<int>(spectrum);
79 const auto x =
m_inWS->points(spectrum);
80 const auto &
y =
m_inWS->y(spectrum).rawData();
81 const auto &e =
m_inWS->e(spectrum).rawData();
84 const auto maxIt = std::max_element(
y.begin(),
y.end());
85 const double height = *maxIt;
86 size_t maxIndex =
static_cast<size_t>(std::distance(
y.begin(), maxIt));
91 size_t leftHalf = maxIndex, rightHalf =
x.size() - maxIndex - 1;
95 for (
auto it = maxIt; it !=
y.end(); ++it) {
97 rightHalf = it - maxIt - 1;
104 for (
auto it = maxIt; it !=
y.begin(); --it) {
106 leftHalf = maxIt - it - 1;
110 g_log.
debug() <<
"Peak in spectrum #" << spectrum <<
" has last bins above 0.5*max at " << leftHalf <<
"\t"
111 << rightHalf <<
"\n";
115 if (rightHalf + leftHalf >= 2) {
118 double fwhm =
x[maxIndex + rightHalf] -
x[maxIndex - leftHalf];
119 double sigma = fwhm / (2. * sqrt(2. * log(2.)));
120 double center =
x[maxIndex];
121 double start = center - 3. * fwhm;
122 double end = center + 3. * fwhm;
124 std::stringstream function;
125 function <<
"name=Gaussian,PeakCentre=";
126 function << center <<
",Height=" <<
height <<
",Sigma=" <<
sigma;
128 g_log.
debug() <<
"Fitting spectrum #" << spectrum <<
" with: " << function.str() <<
"\n";
131 fitAlg->setProperty(
"Function", function.str());
132 fitAlg->setProperty(
"InputWorkspace",
m_inWS);
133 fitAlg->setProperty(
"WorkspaceIndex",
static_cast<int>(spectrum));
134 fitAlg->setProperty(
"StartX", start);
135 fitAlg->setProperty(
"EndX", end);
136 fitAlg->setProperty(
"CreateOutput",
true);
137 fitAlg->setProperty(
"OutputParametersOnly",
true);
138 fitAlg->executeAsChildAlg();
140 const std::string status = fitAlg->getProperty(
"OutputStatus");
143 if (status ==
"success") {
144 m_outWS->cell<
double>(spectrum, 1) = fitResult->cell<
double>(1, 1);
145 m_outWS->cell<
double>(spectrum, 2) = fitResult->cell<
double>(1, 2);
146 m_outWS->cell<
double>(spectrum, 3) = fitResult->cell<
double>(2, 1);
147 m_outWS->cell<
double>(spectrum, 4) = fitResult->cell<
double>(2, 2);
148 m_outWS->cell<
double>(spectrum, 5) = fitResult->cell<
double>(0, 1);
149 m_outWS->cell<
double>(spectrum, 6) = fitResult->cell<
double>(0, 2);
150 m_outWS->cell<
double>(spectrum, 7) = fitResult->cell<
double>(3, 1);
151 m_outWS->cell<std::string>(spectrum, 8) = status;
153 g_log.
debug() <<
"Fit failed in spectrum #" << spectrum <<
". \nReason :" << status
154 <<
". \nSetting the maximum.\n";
155 m_outWS->cell<std::string>(spectrum, 8) =
"fitFailed";
156 m_outWS->cell<
double>(spectrum, 1) =
x[maxIndex];
157 m_outWS->cell<
double>(spectrum, 2) = 0.;
159 m_outWS->cell<
double>(spectrum, 6) = e[maxIndex];
163 g_log.
information() <<
"Found <=3 bins above half maximum in spectrum #" <<
index <<
". Not fitting.\n";
164 m_outWS->cell<std::string>(spectrum, 8) =
"narrowPeak";
165 m_outWS->cell<
double>(spectrum, 1) =
x[maxIndex];
166 m_outWS->cell<
double>(spectrum, 2) = 0.;
168 m_outWS->cell<
double>(spectrum, 6) = e[maxIndex];
171 g_log.
notice() <<
"Negative maximum in spectrum #" << spectrum <<
". Skipping.\n";
172 m_outWS->cell<std::string>(spectrum, 8) =
"negativeMaximum";
182 m_outWS = std::make_shared<TableWorkspace>();
184 const std::vector<std::string> columns = {
"PeakCentre",
"PeakCentreError",
"Sigma",
"SigmaError",
185 "Height",
"HeightError",
"chiSq"};
187 m_outWS->addColumn(
"int",
"WorkspaceIndex");
188 m_outWS->getColumn(0)->setPlotType(1);
189 for (
const auto &column : columns) {
190 m_outWS->addColumn(
"double", column);
192 m_outWS->addColumn(
"str",
"FitStatus");
194 const size_t numberSpectra =
m_inWS->getNumberHistograms();
195 m_progress = std::make_unique<Progress>(
this, 0.0, 1.0, numberSpectra);
197 m_outWS->setRowCount(numberSpectra);
#define DECLARE_ALGORITHM(classname)
std::map< DeltaEMode::Type, std::string > index
#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...
#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_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.
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 property class for workspaces.
Performs Gaussian fits over each spectrum to find the Elastic Peak Position (EPP).
int version() const override
Algorithm's version for identification.
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
void fitGaussian(int64_t)
const std::string category() const override
Algorithm's category for identification.
Mantid::API::ITableWorkspace_sptr m_outWS
Mantid::API::MatrixWorkspace_sptr m_inWS
void init() override
Initialize the algorithm's properties.
void exec() override
Execute the algorithm.
void initWorkspace()
Initializes the output workspace.
std::unique_ptr< Mantid::API::Progress > m_progress
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.
void notice(const std::string &msg)
Logs at notice level.
void information(const std::string &msg)
Logs at information level.
std::shared_ptr< ITableWorkspace > ITableWorkspace_sptr
shared pointer to Mantid::API::ITableWorkspace
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.