Mantid
Loading...
Searching...
No Matches
IndexSXPeaks.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 +
7//----------------------------------------------------------------------
8// Includes
9//----------------------------------------------------------------------
16
17namespace Mantid::Crystal {
18
19// Register the class into the algorithm factory
20DECLARE_ALGORITHM(IndexSXPeaks)
21
22using namespace Geometry;
23using namespace API;
24using namespace Kernel;
25
30 auto mustBePositive = std::make_shared<BoundedValidator<double>>();
31 mustBePositive->setLower(0.0);
32
33 auto reasonable_angle = std::make_shared<BoundedValidator<double>>();
34 reasonable_angle->setLower(5.0);
35 reasonable_angle->setUpper(175.0);
36
38 std::make_unique<WorkspaceProperty<Mantid::DataObjects::PeaksWorkspace>>("PeaksWorkspace", "", Direction::InOut),
39 "Input Peaks Workspace");
40
41 declareProperty(std::make_unique<PropertyWithValue<double>>("a", -1.0, mustBePositive->clone(), Direction::Input),
42 "Lattice parameter a");
43
44 declareProperty(std::make_unique<PropertyWithValue<double>>("b", -1.0, mustBePositive->clone(), Direction::Input),
45 "Lattice parameter b");
46
47 declareProperty(std::make_unique<PropertyWithValue<double>>("c", -1.0, std::move(mustBePositive), Direction::Input),
48 "Lattice parameter c");
49
51 std::make_unique<PropertyWithValue<double>>("alpha", -1.0, reasonable_angle->clone(), Direction::Input),
52 "Lattice parameter alpha");
53
55 std::make_unique<PropertyWithValue<double>>("beta", -1.0, reasonable_angle->clone(), Direction::Input),
56 "Lattice parameter beta");
57
59 std::make_unique<PropertyWithValue<double>>("gamma", -1.0, std::move(reasonable_angle), Direction::Input),
60 "Lattice parameter gamma");
61
62 declareProperty(std::make_unique<ArrayProperty<int>>("PeakIndices"),
63 "Index of the peaks in the table workspace to be used. If no "
64 "index are provided, all will be used.");
65
66 declareProperty("dTolerance", 0.01, "Tolerance for peak positions in d-spacing");
67
68 std::vector<int> extents(6, 0);
69 const int range = 20;
70 extents[0] = -range;
71 extents[1] = range;
72 extents[2] = -range;
73 extents[3] = range;
74 extents[4] = -range;
75 extents[5] = range;
76 declareProperty(std::make_unique<ArrayProperty<int>>("SearchExtents", std::move(extents)),
77 "A comma separated list of min, max for each of H, K and L,\n"
78 "Specifies the search extents applied for H K L values "
79 "associated with the peaks.");
80}
81
88void IndexSXPeaks::cullHKLs(std::vector<PeakCandidate> &peakCandidates, Mantid::Geometry::UnitCell &unitcell) {
89 size_t npeaks = peakCandidates.size();
90 for (std::size_t p = 0; p < npeaks; p++) {
91 for (std::size_t q = 0; q < npeaks; q++) {
92 if (p == q) // Don't do a self comparison
93 {
94 continue;
95 }
96 peakCandidates[p].clean(peakCandidates[q], unitcell,
97 0.5 * M_PI / 180.0); // Half a degree tolerance
98 }
99 }
100}
101
107void IndexSXPeaks::validateNotColinear(const std::vector<PeakCandidate> &peakCandidates) const {
108 // Find two non-colinear peaks
109 bool all_collinear = true;
110 size_t npeaks = peakCandidates.size();
111 for (std::size_t i = 0; i < npeaks; i++) {
112 for (std::size_t j = i; j < npeaks; j++) {
113 double anglerad = peakCandidates[i].angle(peakCandidates[j]);
114 if (anglerad > 2.0 * M_PI / 180.0 && anglerad < 178.0 * M_PI / 180.0) {
115 all_collinear = false;
116 break;
117 }
118 }
119 }
120 // Throw if all collinear
121 if (all_collinear) {
122 throw std::runtime_error("Angles between all pairs of peaks are too small");
123 }
124}
125
131 using namespace Mantid::DataObjects;
132 std::vector<int> peakindices = getProperty("PeakIndices");
133
134 PeaksWorkspace_sptr ws = this->getProperty("PeaksWorkspace");
135
136 // Need a least two peaks
137 std::size_t npeaks = peakindices.size();
138 if (npeaks > size_t(ws->getNumberPeaks())) {
139 throw std::runtime_error("Cannot have more peaks indices than actual peaks");
140 }
141 if (npeaks == 1 || ws->getNumberPeaks() < 2) {
142 throw std::runtime_error("At least 2 peaks are required for this algorithm to work");
143 }
144 if (npeaks == 0) {
145 // If the user provides no peaks we default to use all the available peaks.
146 npeaks = ws->getNumberPeaks();
147 peakindices.reserve(npeaks);
148 for (int i = 1; i <= int(npeaks); i++) // create indexes corresponding to all peak indexes
149 {
150 peakindices.emplace_back(i);
151 }
152 g_log.information("No peak indexes provided. Algorithm will use all peaks "
153 "in the workspace for the calculation.");
154 }
155
156 // Get the effective unit cell
157 double a = getProperty("a");
158 double b = getProperty("b");
159 double c = getProperty("c");
160 double alpha = getProperty("alpha");
161 double beta = getProperty("beta");
162 double gamma = getProperty("gamma");
163
164 std::vector<int> extents = getProperty("SearchExtents");
165 if (extents.size() != 6) {
166 std::stringstream stream;
167 stream << "Expected 6 elements for the extents. Got: " << extents.size();
168 throw std::runtime_error(stream.str());
169 }
170
171 // Create the Unit-Cell.
172 Mantid::Geometry::UnitCell unitcell(a, b, c, alpha, beta, gamma);
173
174 std::vector<PeakCandidate> peaks;
175
176 for (std::size_t i = 0; i < npeaks; i++) {
177 int row = peakindices[i] - 1;
178 if (row < 0) {
179 throw std::runtime_error("Cannot have a peak index < 0.");
180 }
181 IPeak &peak = ws->getPeak(row);
182 V3D Qs = peak.getQSampleFrame() / (2.0 * M_PI);
183 peaks.emplace_back(Qs[0], Qs[1], Qs[2]);
184 }
185 assert(peaks.size() >= 2);
186
187 // Sanity check the generated peaks.
188 validateNotColinear(peaks);
189
190 // Generate HKL possibilities for each peak.
191 double dtol = getProperty("dTolerance");
192 Progress prog(this, 0.0, 1.0, 4);
193 for (int h = extents[0]; h < extents[1]; h++) {
194 for (int k = extents[2]; k < extents[3]; k++) {
195 for (int l = extents[4]; l < extents[5]; l++) {
196 double dspacing = unitcell.d(h, k, l); // Create a fictional d spacing
197 for (std::size_t p = 0; p < npeaks; p++) {
198 double dSpacingPeaks = peaks[p].getdSpacing();
199 if (std::abs(dspacing - dSpacingPeaks) < dtol)
200 peaks[p].addHKL(h, k, l); // If the peak position and the fictional
201 // d spacing are within tolerance, add it
202 }
203 }
204 }
205 }
206 prog.report(); // 1st Progress report.
207
208 cullHKLs(peaks, unitcell);
209
210 prog.report(); // 2nd progress report.
211 peaks[0].setFirst(); // On the first peak, now only the first candidate hkl is
212 // considered, others are erased,
213 // This means the design space of possible peak-hkl alignments has been
214 // reduced, will improve future refinements.
215
216 cullHKLs(peaks, unitcell);
217 prog.report(); // 3rd progress report.
218
219 peaks[1].setFirst();
220
221 cullHKLs(peaks, unitcell);
222 prog.report(); // 4th progress report.
223
224 // Now we can index the input/output peaks workspace
225 // If there are peak indexes uses those to find actual peaks in the workspace
226 // and overrite HKL
227 for (std::size_t i = 0; i < npeaks; i++) {
228 int row = 0;
229 try {
230 row = peakindices[i] - 1;
231 IPeak &peak = ws->getPeak(row);
232 const V3D hkl = peaks[i].getHKL();
233 peak.setHKL(hkl);
234 std::stringstream stream;
235 stream << "Peak Index: " << row << " HKL: " << hkl;
236 g_log.information(stream.str());
237 } catch (std::logic_error &) {
238 std::stringstream msg;
239 msg << "Peak Index: " << row + 1 << " cannot be assigned a single HKL set.";
240 g_log.warning(msg.str());
241 continue;
242 }
243 }
244}
245} // namespace Mantid::Crystal
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
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
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A property class for workspaces.
void validateNotColinear(const std::vector< PeakCandidate > &peakCandidates) const
Check that not all peaks are colinear and throw if they are not.
void init() override
Initialisation method.
void exec() override
Executes the algorithm.
void cullHKLs(std::vector< PeakCandidate > &peakCandidates, Mantid::Geometry::UnitCell &unitcell)
Culling method to direct the removal of hkl values off peaks where they cannot sit.
Structure describing a single-crystal peak.
Definition: IPeak.h:26
virtual void setHKL(double H, double K, double L)=0
virtual Mantid::Kernel::V3D getQSampleFrame() const =0
Class to implement unit cell of crystals.
Definition: UnitCell.h:44
double d(double h, double k, double l) const
Return d-spacing ( ) for a given h,k,l coordinate.
Definition: UnitCell.cpp:700
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
The concrete, templated class for properties.
Class for 3D vectors.
Definition: V3D.h:34
std::shared_ptr< PeaksWorkspace > PeaksWorkspace_sptr
Typedef for a shared pointer to a peaks workspace.
static constexpr double h
Planck constant in J*s.
@ InOut
Both an input & output workspace.
Definition: Property.h:55
@ Input
An input workspace.
Definition: Property.h:53