-
Notifications
You must be signed in to change notification settings - Fork 122
/
FindUBUsingIndexedPeaks.cpp
236 lines (204 loc) · 8.54 KB
/
FindUBUsingIndexedPeaks.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
// NScD Oak Ridge National Laboratory, European Spallation Source,
// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidCrystal/FindUBUsingIndexedPeaks.h"
#include "MantidAPI/IPeaksWorkspace.h"
#include "MantidAPI/Sample.h"
#include "MantidDataObjects/LeanElasticPeaksWorkspace.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidGeometry/Crystal/IndexingUtils.h"
#include "MantidGeometry/Crystal/OrientedLattice.h"
#include "MantidKernel/BoundedValidator.h"
namespace Mantid::Crystal {
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(FindUBUsingIndexedPeaks)
using namespace Mantid::Kernel;
using namespace Mantid::API;
using namespace Mantid::DataObjects;
using namespace Mantid::Geometry;
const int MIN_INDEXED_PEAKS = 3;
/** Initialize the algorithm's properties.
*/
void FindUBUsingIndexedPeaks::init() {
declareProperty(std::make_unique<WorkspaceProperty<IPeaksWorkspace>>("PeaksWorkspace", "", Direction::InOut),
"Input Peaks Workspace");
auto mustBePositive = std::make_shared<BoundedValidator<double>>();
mustBePositive->setLower(0.0);
declareProperty("Tolerance", 0.1, mustBePositive, "Indexing Tolerance (0.1)");
declareProperty("ToleranceForSatellite", 0.1, mustBePositive, "Indexing Tolerance for satellite (0.1)");
declareProperty("CommonUBForAll", false, "Used when evaluating the uncertainty of modHKL");
}
/** Execute the algorithm.
*/
void FindUBUsingIndexedPeaks::exec() {
IPeaksWorkspace_sptr ws = getProperty("PeaksWorkspace");
const int n_peaks = ws->getNumberPeaks();
std::vector<V3D> q_vectors;
std::vector<V3D> hkl_vectors;
std::vector<V3D> mnp_vectors;
double tolerance = getProperty("Tolerance");
bool commonUB = getProperty("CommonUBForAll");
int ModDim = 0;
int MaxOrder = 0;
bool CrossTerm = false;
DblMatrix errorHKL(3, 3);
q_vectors.reserve(n_peaks);
hkl_vectors.reserve(n_peaks);
mnp_vectors.reserve(n_peaks);
size_t indexed_count = 0;
std::unordered_set<int> run_numbers;
for (int i = 0; i < n_peaks; i++) {
const IPeak &peak = ws->getPeak(i);
run_numbers.insert(peak.getRunNumber());
V3D hkl(peak.getIntHKL());
V3D mnp(peak.getIntMNP());
int maxOrder = mnp.maxCoeff();
if (maxOrder > MaxOrder)
MaxOrder = maxOrder;
if (mnp[0] != 0 && ModDim == 0)
ModDim = 1;
if (mnp[1] != 0 && ModDim == 1)
ModDim = 2;
if (mnp[2] != 0 && ModDim == 2)
ModDim = 3;
if (mnp[0] * mnp[1] != 0 || mnp[1] * mnp[2] != 0 || mnp[2] * mnp[0] != 0)
CrossTerm = true;
if (isPeakIndexed(peak)) {
q_vectors.emplace_back(peak.getQSampleFrame());
hkl_vectors.emplace_back(hkl);
mnp_vectors.emplace_back(mnp);
indexed_count++;
}
}
// too few indexed peaks to work with
if (indexed_count < MIN_INDEXED_PEAKS) {
throw std::runtime_error("At least three linearly independent indexed peaks are needed.");
}
Matrix<double> UB(3, 3, false);
Matrix<double> modUB(3, 3, false);
std::vector<double> sigabc(7);
std::vector<double> sigq(3);
IndexingUtils::Optimize_6dUB(UB, modUB, hkl_vectors, mnp_vectors, ModDim, q_vectors, sigabc, sigq);
if (!IndexingUtils::CheckUB(UB)) // UB not found correctly
{
g_log.notice(std::string("Found Invalid UB...peaks used might not be linearly independent"));
g_log.notice(std::string("UB NOT SAVED."));
} else // tell user how many would be indexed
{ // from the full list of peaks, and
q_vectors.clear(); // save the UB in the sample
q_vectors.reserve(n_peaks);
for (int i = 0; i < n_peaks; i++) {
q_vectors.emplace_back(ws->getPeak(i).getQSampleFrame());
}
int num_indexed = IndexingUtils::NumberIndexed(UB, q_vectors, tolerance);
int sate_indexed = 0;
double satetolerance = getProperty("ToleranceForSatellite");
if (ModDim > 0) {
for (int run : run_numbers) {
std::vector<V3D> run_q_vectors;
std::vector<V3D> run_fhkl_vectors;
std::vector<V3D> run_hkl_vectors;
std::vector<V3D> run_mnp_vectors;
size_t run_indexed = 0;
for (int i = 0; i < n_peaks; i++) {
const IPeak &peak = ws->getPeak(i);
if (peak.getRunNumber() == run) {
if (isPeakIndexed(peak))
run_indexed++;
}
}
g_log.notice() << "Number of Indexed Peaks in Run " << run << " is " << run_indexed << "\n";
run_q_vectors.reserve(run_indexed);
run_hkl_vectors.reserve(run_indexed);
run_mnp_vectors.reserve(run_indexed);
run_fhkl_vectors.reserve(run_indexed);
for (int i = 0; i < n_peaks; i++) {
const IPeak &peak = ws->getPeak(i);
if (peak.getRunNumber() == run) {
V3D hkl(peak.getIntHKL());
V3D mnp(peak.getIntMNP());
if (isPeakIndexed(peak)) {
run_q_vectors.emplace_back(peak.getQSampleFrame());
run_hkl_vectors.emplace_back(hkl);
run_mnp_vectors.emplace_back(mnp);
}
}
}
if (run_indexed < 3)
continue;
OrientedLattice run_lattice;
if (commonUB) {
run_lattice.setUB(UB);
run_lattice.setModUB(modUB);
run_lattice.setError(sigabc[0], sigabc[1], sigabc[2], sigabc[3], sigabc[4], sigabc[5]);
} else {
Matrix<double> run_UB(3, 3, false);
Matrix<double> run_modUB(3, 3, false);
std::vector<double> run_sigabc(7);
std::vector<double> run_sigq(3);
IndexingUtils::Optimize_6dUB(run_UB, run_modUB, run_hkl_vectors, run_mnp_vectors, ModDim, run_q_vectors,
run_sigabc, run_sigq);
run_lattice.setUB(run_UB);
run_lattice.setModUB(run_modUB);
run_lattice.setError(run_sigabc[0], run_sigabc[1], run_sigabc[2], run_sigabc[3], run_sigabc[4],
run_sigabc[5]);
}
g_log.notice() << run_lattice << "\n";
double average_error = 0.;
IndexingUtils::CalculateMillerIndices(run_lattice.getUB(), run_q_vectors, 1.0, run_fhkl_vectors, average_error);
for (size_t i = 0; i < run_indexed; i++) {
if (IndexingUtils::ValidIndex(run_fhkl_vectors[i], tolerance))
continue;
V3D fhkl(run_fhkl_vectors[i]);
for (int j = 0; j < 3; j++) {
if (run_mnp_vectors[i][j] != 0) {
fhkl -= run_lattice.getModVec(j) * run_mnp_vectors[i][j];
if (IndexingUtils::ValidIndex(fhkl, satetolerance)) {
sate_indexed++;
V3D errhkl = fhkl - run_hkl_vectors[i];
errhkl = errhkl.absoluteValue();
for (int k = 0; k < 3; k++)
errorHKL[k][j] += errhkl[k];
}
}
}
}
}
}
std::stringstream stream;
stream << "New UB will index " << num_indexed << " main Peaks with tolerance " << tolerance << " and "
<< sate_indexed << " Satellite Peaks with tolerance " << satetolerance << " ,out of " << n_peaks
<< " Peaks \n";
g_log.notice(stream.str());
auto o_lattice = std::make_unique<OrientedLattice>();
o_lattice->setUB(UB);
o_lattice->setModUB(modUB);
o_lattice->setError(sigabc[0], sigabc[1], sigabc[2], sigabc[3], sigabc[4], sigabc[5]);
double ind_count_inv = 1. / static_cast<double>(indexed_count);
errorHKL *= ind_count_inv;
o_lattice->setErrorModHKL(errorHKL);
o_lattice->setMaxOrder(MaxOrder);
o_lattice->setCrossTerm(CrossTerm);
// Show the modified lattice parameters
logLattice(*o_lattice, ModDim);
ws->mutableSample().setOrientedLattice(std::move(o_lattice));
}
}
void FindUBUsingIndexedPeaks::logLattice(const OrientedLattice &o_lattice, const int &ModDim) {
// Show the modified lattice parameters
g_log.notice() << o_lattice << "\n";
g_log.notice() << "Modulation Dimension is: " << ModDim << "\n";
for (int i = 0; i < ModDim; i++) {
g_log.notice() << "Modulation Vector " << i + 1 << ": " << o_lattice.getModVec(i) << "\n";
g_log.notice() << "Modulation Vector " << i + 1 << " error: " << o_lattice.getVecErr(i) << "\n";
}
}
bool FindUBUsingIndexedPeaks::isPeakIndexed(const IPeak &peak) {
const V3D hkl(peak.getIntHKL());
const V3D mnp(peak.getIntMNP());
return (IndexingUtils::ValidIndex(hkl, 1.0) || IndexingUtils::ValidIndex(mnp, 1.0));
}
} // namespace Mantid::Crystal