Mantid
Loading...
Searching...
No Matches
GeneralisedSecondDifference.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
11#include "MantidIndexing/Extract.h"
12#include "MantidIndexing/IndexInfo.h"
15
16#include <numeric>
17#include <sstream>
18
19namespace Mantid::Algorithms {
20
21// Register the class into the algorithm factory
22DECLARE_ALGORITHM(GeneralisedSecondDifference)
23
24using namespace Kernel;
25using namespace API;
26
29
30 // Input and output workspaces
31 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("InputWorkspace", "", Direction::Input),
32 "Name of the input workspace");
33 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("OutputWorkspace", "", Direction::Output),
34 "The name of the workspace to be created as the output of the algorithm");
35
36 auto mustBePositive = std::make_shared<BoundedValidator<int>>();
37 mustBePositive->setLower(0);
38 declareProperty("M", 0, mustBePositive,
39 "The number of points for averaging, i.e. summing will be done in the\n"
40 "range [y(i-m),y(i+m)]");
41 declareProperty("Z", 0, mustBePositive, "The number of iteration steps in the averaging procedure");
42 declareProperty("WorkspaceIndexMin", 0, mustBePositive, "Lower bound of the spectrum range (default 0)");
43 declareProperty("WorkspaceIndexMax", 0, mustBePositive, "Upper bound of the spectrum range (default workspace max)");
44}
45
51 // Message stream used for logger
52 std::ostringstream message;
53
54 // Get some properties
55 MatrixWorkspace_const_sptr inputWS = getProperty("InputWorkspace");
56 int spec_min = getProperty("WorkspaceIndexMin");
57 int spec_max = getProperty("WorkspaceIndexMax");
58 auto n_hists = static_cast<int>(inputWS->getNumberHistograms());
59
60 if (spec_min == 0 && spec_max == 0) // Values per default, take all spectra
61 spec_max = n_hists - 1;
62
63 if (spec_min > spec_max)
64 std::swap(spec_min, spec_max);
65
66 if (spec_max > n_hists) {
67 message << "WorkspaceIndexMax " << spec_max << " > number of histograms, WorkspaceIndexMax reset to "
68 << (n_hists - 1);
69 g_log.information(message.str());
70 message.str("");
71 spec_max = n_hists - 1;
72 }
73
74 // Get some more input fields
75 m_z = getProperty("Z");
76 m_m = getProperty("M");
77 const int n_av = m_z * m_m + 1;
78
79 // Calculate the Cij and Cij^2 coefficients
81
82 const int n_points = static_cast<int>(inputWS->y(0).size()) - 2 * n_av;
83 if (n_points < 1) {
84 throw std::invalid_argument("Invalid (M,Z) values");
85 }
86 // Create OuputWorkspace
87 auto out = DataObjects::create<HistoWorkspace>(*inputWS, Indexing::extract(inputWS->indexInfo(), spec_min, spec_max),
88 HistogramData::BinEdges(n_points + 1));
89
90 const int nsteps = 2 * n_av + 1;
91
92 std::shared_ptr<API::Progress> progress = std::make_shared<API::Progress>(this, 0.0, 1.0, (spec_max - spec_min));
93 for (int i = spec_min; i <= spec_max; i++) {
94 int out_index = i - spec_min;
95 const auto &refX = inputWS->x(i);
96 const auto &refY = inputWS->y(i);
97 const auto &refE = inputWS->e(i);
98 auto &outX = out->mutableX(out_index);
99 auto &outY = out->mutableY(out_index);
100 auto &outE = out->mutableE(out_index);
101
102 std::copy(refX.begin() + n_av, refX.end() - n_av, outX.begin());
103
104 auto itInY = refY.cbegin();
105 auto itOutY = outY.begin();
106 auto itInE = refE.cbegin();
107 auto itOutE = outE.begin();
108 for (; itOutY != outY.end(); ++itOutY, ++itInY, ++itOutE, ++itInE) {
109 // Calculate \sum_{j}Cij.Y(j)
110 (*itOutY) = std::inner_product(itInY, itInY + nsteps, m_Cij.begin(), 0.0);
111 // Calculate the error bars sqrt(\sum_{j}Cij^2.E^2)
112 double err2 = std::inner_product(itInE, itInE + nsteps, m_Cij2.cbegin(), 0.0);
113 (*itOutE) = sqrt(err2);
114 }
115 progress->report();
116 }
117 setProperty("OutputWorkspace", std::move(out));
118}
123 int zz = 0;
124 int max_index_prev = 1;
125 int n_el_prev = 3;
126 std::vector<double> previous(n_el_prev);
127 previous[0] = 1;
128 previous[1] = -2;
129 previous[2] = 1;
130
131 if (m_z == 0) //
132 {
133 m_Cij.resize(3);
134 m_Cij.assign(previous.begin(), previous.end());
135 m_Cij2.resize(3);
136 std::transform(m_Cij.cbegin(), m_Cij.cend(), m_Cij2.begin(), VectorHelper::Squares<double>());
137 return;
138 }
139 std::vector<double> next;
140 // Calculate the Cij iteratively.
141 do {
142 zz++;
143 int max_index = zz * m_m + 1;
144 int n_el = 2 * max_index + 1;
145 next.resize(n_el);
146 std::fill(next.begin(), next.end(), 0.0);
147 for (int i = 0; i < n_el; ++i) {
148 int delta = -max_index + i;
149 for (int l = delta - m_m; l <= delta + m_m; l++) {
150 int index = l + max_index_prev;
151 if (index >= 0 && index < n_el_prev)
152 next[i] += previous[index];
153 }
154 }
155 previous = next;
156 max_index_prev = max_index;
157 n_el_prev = n_el;
158 } while (zz != m_z);
159
160 m_Cij.resize(2 * m_z * m_m + 3);
161 m_Cij.assign(previous.begin(), previous.end());
162 m_Cij2.resize(2 * m_z * m_m + 3);
163 std::transform(m_Cij.cbegin(), m_Cij.cend(), m_Cij2.begin(), VectorHelper::Squares<double>());
164}
165
166} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
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
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
Definition: Algorithm.cpp:231
A property class for workspaces.
int m_z
Contains the value of the property z.
std::vector< double > m_Cij2
Vector that contains the prefactor coefficients Cij^2 in the range [-zm-1,zm+1].
std::vector< double > m_Cij
Vector that contains the prefactor coefficients Cij in the range [-zm-1,zm+1].
int m_m
Contains the value of the property m.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54