Mantid
Loading...
Searching...
No Matches
FindEPP.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 +
10
11#include <cmath>
12#include <sstream>
13
14namespace Mantid::Algorithms {
15
16using namespace Mantid::Kernel;
17using namespace Mantid::API;
18using namespace Mantid::DataObjects;
19
20// Register the algorithm into the AlgorithmFactory
21DECLARE_ALGORITHM(FindEPP)
22
23//----------------------------------------------------------------------------------------------
24
25
26const std::string FindEPP::name() const { return "FindEPP"; }
27
29int FindEPP::version() const { return 2; }
30
32const std::string FindEPP::category() const { return "Workflow\\MLZ\\TOFTOF;Utility"; }
33
35const std::string FindEPP::summary() const {
36 return "Performs Gaussian fits over each spectrum to find the Elastic Peak "
37 "Position (EPP).";
38}
39
40//----------------------------------------------------------------------------------------------
44 declareProperty(std::make_unique<WorkspaceProperty<API::MatrixWorkspace>>("InputWorkspace", "", Direction::Input),
45 "An input workspace.");
46 declareProperty(std::make_unique<WorkspaceProperty<API::ITableWorkspace>>("OutputWorkspace", "", Direction::Output),
47 "An output workspace.");
48}
49
50//----------------------------------------------------------------------------------------------
54 m_inWS = getProperty("InputWorkspace");
55
57
58 auto numberspectra = static_cast<int64_t>(m_inWS->getNumberHistograms());
59
60 // Loop over spectra
62 for (int64_t index = 0; index < numberspectra; ++index) {
66 }
68
69 setProperty("OutputWorkspace", m_outWS);
70}
71
72/* Call Fit as child algorithm for each spectra
73 * @param index : the workspace index
74 */
76 auto spectrum = static_cast<size_t>(index);
77 m_outWS->cell<int>(spectrum, 0) = static_cast<int>(spectrum);
78
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();
82
83 // Find the maximum value and it's index
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));
87
88 if (height > 0) {
89 // Find how many bins are around maximum, that are above half-maximum
90 // Initialize the distances of the half-maxima bins from maximum
91 size_t leftHalf = maxIndex, rightHalf = x.size() - maxIndex - 1;
92
93 // Find the first bin on the right side of maximum, that drops below
94 // half-maximum
95 for (auto it = maxIt; it != y.end(); ++it) {
96 if (*it < 0.5 * height) {
97 rightHalf = it - maxIt - 1;
98 break;
99 }
100 }
101
102 // Find the first bin on the left side of maximum, that drops below
103 // half-maximum
104 for (auto it = maxIt; it != y.begin(); --it) {
105 if (*it < 0.5 * height) {
106 leftHalf = maxIt - it - 1;
107 break;
108 }
109 }
110 g_log.debug() << "Peak in spectrum #" << spectrum << " has last bins above 0.5*max at " << leftHalf << "\t"
111 << rightHalf << "\n";
112
113 // We want to fit only if there are at least 3 bins (including the maximum
114 // itself) above half-maximum
115 if (rightHalf + leftHalf >= 2) {
116
117 // Prepare the initial parameters for the fit
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;
123
124 std::stringstream function;
125 function << "name=Gaussian,PeakCentre=";
126 function << center << ",Height=" << height << ",Sigma=" << sigma;
127
128 g_log.debug() << "Fitting spectrum #" << spectrum << " with: " << function.str() << "\n";
129
130 auto fitAlg = createChildAlgorithm("Fit", 0., 0., false);
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();
139
140 const std::string status = fitAlg->getProperty("OutputStatus");
141 ITableWorkspace_sptr fitResult = fitAlg->getProperty("OutputParameters");
142
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;
152 } else {
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.;
158 m_outWS->cell<double>(spectrum, 5) = height;
159 m_outWS->cell<double>(spectrum, 6) = e[maxIndex];
160 }
161
162 } else {
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.;
167 m_outWS->cell<double>(spectrum, 5) = height;
168 m_outWS->cell<double>(spectrum, 6) = e[maxIndex];
169 }
170 } else {
171 g_log.notice() << "Negative maximum in spectrum #" << spectrum << ". Skipping.\n";
172 m_outWS->cell<std::string>(spectrum, 8) = "negativeMaximum";
173 }
174 m_progress->report();
175}
176
181
182 m_outWS = std::make_shared<TableWorkspace>();
183
184 const std::vector<std::string> columns = {"PeakCentre", "PeakCentreError", "Sigma", "SigmaError",
185 "Height", "HeightError", "chiSq"};
186
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);
191 }
192 m_outWS->addColumn("str", "FitStatus");
193
194 const size_t numberSpectra = m_inWS->getNumberHistograms();
195 m_progress = std::make_unique<Progress>(this, 0.0, 1.0, numberSpectra);
196
197 m_outWS->setRowCount(numberSpectra);
198}
199
200} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double sigma
Definition: GetAllEi.cpp:156
double height
Definition: GetAllEi.cpp:155
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
#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_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
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 property class for workspaces.
Performs Gaussian fits over each spectrum to find the Elastic Peak Position (EPP).
Definition: FindEPP.h:20
int version() const override
Algorithm's version for identification.
Definition: FindEPP.cpp:29
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
Definition: FindEPP.cpp:35
void fitGaussian(int64_t)
Definition: FindEPP.cpp:75
const std::string category() const override
Algorithm's category for identification.
Definition: FindEPP.cpp:32
Mantid::API::ITableWorkspace_sptr m_outWS
Definition: FindEPP.h:34
Mantid::API::MatrixWorkspace_sptr m_inWS
Definition: FindEPP.h:33
void init() override
Initialize the algorithm's properties.
Definition: FindEPP.cpp:43
void exec() override
Execute the algorithm.
Definition: FindEPP.cpp:53
void initWorkspace()
Initializes the output workspace.
Definition: FindEPP.cpp:180
std::unique_ptr< Mantid::API::Progress > m_progress
Definition: FindEPP.h:35
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
void notice(const std::string &msg)
Logs at notice level.
Definition: Logger.cpp:95
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
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.
Definition: MultiThreaded.h:22
STL namespace.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54