Mantid
Loading...
Searching...
No Matches
IndexPeaks.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 +
9#include "MantidAPI/Sample.h"
18
19#include <boost/optional/optional.hpp>
20
21#include <algorithm>
22
23namespace Mantid::Crystal {
24// Register the algorithm into the AlgorithmFactory
25DECLARE_ALGORITHM(IndexPeaks)
26
27using API::IPeaksWorkspace;
29using Geometry::IndexingUtils;
30using Geometry::IPeak;
31using Geometry::OrientedLattice;
33using Kernel::Logger;
34using Kernel::V3D;
35
36namespace {
37const auto OPTIMIZE_UB_ATTEMPTS{4};
38
39namespace Prop {
40const std::string PEAKSWORKSPACE{"PeaksWorkspace"};
41const std::string TOLERANCE{"Tolerance"};
42const std::string SATE_TOLERANCE{"ToleranceForSatellite"};
43const std::string ROUNDHKLS{"RoundHKLs"};
44const std::string COMMONUB{"CommonUBForAll"};
45const std::string SAVEMODINFO{"SaveModulationInfo"};
46const std::string AVERAGE_ERR{"AverageError"};
47const std::string NUM_INDEXED{"NumIndexed"};
48const std::string MAIN_NUM_INDEXED{"MainNumIndexed"};
49const std::string SATE_NUM_INDEXED{"SateNumIndexed"};
50const std::string MAIN_ERR{"MainError"};
51const std::string SATE_ERR{"SatelliteError"};
52
53struct SatelliteIndexingArgs {
54 const double tolerance;
55 const int maxOrder;
56 const std::vector<V3D> modVectors;
57 const bool crossTerms;
58};
59
60struct IndexPeaksArgs {
66 static IndexPeaksArgs parse(const API::Algorithm &alg) {
67 const IPeaksWorkspace_sptr peaksWS = alg.getProperty(PEAKSWORKSPACE);
68 const int maxOrderFromAlg = alg.getProperty(ModulationProperties::MaxOrder);
69
70 // Init variables
71 int maxOrderToUse{0};
72 std::vector<V3D> modVectorsToUse;
73 modVectorsToUse.reserve(3);
74 bool crossTermToUse{false};
75
76 // Parse mod vectors
77 modVectorsToUse = addModulationVectors(alg.getProperty(ModulationProperties::ModVector1),
78 alg.getProperty(ModulationProperties::ModVector2),
79 alg.getProperty(ModulationProperties::ModVector3));
80 // check the 3 mod vectors added from properties
81 modVectorsToUse = validModulationVectors(modVectorsToUse[0], modVectorsToUse[1], modVectorsToUse[2]);
82
83 // deal with case: max order > 0 and no mod vector is specified
84 if (maxOrderFromAlg > 0 && modVectorsToUse.size() == 0) {
85 // Max Order is larger than zero but no modulated vector specified
86 // Assume that the caller method will handle this
87 maxOrderToUse = maxOrderFromAlg;
88 } else if (maxOrderFromAlg == 0 && modVectorsToUse.size() == 0) {
89 // Use lattice definitions if they exist
90 const auto &lattice = peaksWS->sample().getOrientedLattice();
91 crossTermToUse = lattice.getCrossTerm();
92 maxOrderToUse = lattice.getMaxOrder(); // the lattice can return a 0 here
93
94 // if lattice has maxOrder, we will use the modVec from it, otherwise
95 // stick to the input got from previous assignment
96 if (maxOrderToUse > 0) {
97 modVectorsToUse = validModulationVectors(lattice.getModVec(0), lattice.getModVec(1), lattice.getModVec(2));
98 }
99 } else {
100 // Use user specified
101 // default behavior: map everything automatically
102 maxOrderToUse = maxOrderFromAlg;
103 crossTermToUse = alg.getProperty(ModulationProperties::CrossTerms);
104 }
105
106 return {peaksWS,
107 alg.getProperty(TOLERANCE),
108 alg.getProperty(ROUNDHKLS),
109 alg.getProperty(COMMONUB),
110 alg.getProperty(SAVEMODINFO),
111 SatelliteIndexingArgs{alg.getProperty(SATE_TOLERANCE), maxOrderToUse, modVectorsToUse, crossTermToUse}};
112 }
113
115 const double mainTolerance;
116 const bool roundHKLs;
117 const bool commonUB;
119 const SatelliteIndexingArgs satellites;
120};
121} // namespace Prop
122
126struct PeakIndexingStats {
127 PeakIndexingStats &operator+=(const PeakIndexingStats &rhs) {
128 numIndexed += rhs.numIndexed;
129 error += rhs.error;
130 return *this;
131 }
133 double error{0.0};
134};
135
140struct CombinedIndexingStats {
141 CombinedIndexingStats &operator+=(const CombinedIndexingStats &rhs) {
142 main += rhs.main;
143 satellites += rhs.satellites;
144 return *this;
145 }
146
148 inline int totalNumIndexed() const { return main.numIndexed + satellites.numIndexed; }
150 inline double mainError() const {
151 if (main.numIndexed == 0)
152 return 0.0;
153 return main.error / main.numIndexed;
154 }
156 inline double satelliteError() const {
157 if (satellites.numIndexed == 0)
158 return 0.0;
159 return satellites.error / satellites.numIndexed;
160 }
162 inline double averageError() const { return (main.error + satellites.error) / totalNumIndexed(); }
163
164 PeakIndexingStats main;
165 PeakIndexingStats satellites;
166};
167
174DblMatrix optimizeUBMatrix(const DblMatrix &ubOrig, const std::vector<V3D> &qSample, const double tolerance) {
175 DblMatrix optimizedUB(ubOrig);
176
177 double errorAtStart{0.0};
178 std::vector<V3D> millerIndices;
179 millerIndices.reserve(qSample.size());
180 const int numIndexedAtStart =
181 IndexingUtils::CalculateMillerIndices(optimizedUB, qSample, tolerance, millerIndices, errorAtStart);
182
183 if (numIndexedAtStart < 3) {
184 // can't optimize without at least 3 indexed peaks
185 return optimizedUB;
186 }
187
188 for (auto i = 0; i < OPTIMIZE_UB_ATTEMPTS; ++i) {
189 try {
190 // optimization requires rounded indices
191 IndexingUtils::RoundHKLs(millerIndices);
192 IndexingUtils::Optimize_UB(optimizedUB, millerIndices, qSample);
193 } catch (...) {
194 // If there is any problem, such as too few
195 // independent peaks, just use the original UB
196 optimizedUB = ubOrig;
197 break;
198 }
199 double errorInLoop{0.0};
200 const int numIndexedInLoop =
201 IndexingUtils::CalculateMillerIndices(optimizedUB, qSample, tolerance, millerIndices, errorInLoop);
202 if (numIndexedInLoop < numIndexedAtStart) // use the original UB
203 break;
204 }
205 return optimizedUB;
206}
207
209using IndexedSatelliteInfo = std::tuple<V3D, V3D, double>;
210
234boost::optional<IndexedSatelliteInfo> indexSatellite(const V3D &mainHKL, const int maxOrder,
235 const std::vector<V3D> &modVectors, const double tolerance,
236 const bool crossTerms) {
237 const auto offsets = generateOffsetVectors(modVectors, maxOrder, crossTerms);
238 bool foundSatellite{false};
239 V3D indexedIntHKL, indexedMNP;
240 for (const auto &mnpOffset : offsets) {
241 const auto candidateIntHKL = mainHKL - std::get<3>(mnpOffset);
242 const V3D candidateMNP{std::get<0>(mnpOffset), std::get<1>(mnpOffset), std::get<2>(mnpOffset)};
243 if (IndexingUtils::ValidIndex(candidateIntHKL, tolerance)) {
244 indexedIntHKL = candidateIntHKL;
245 indexedMNP = candidateMNP;
246 foundSatellite = true;
247 // we deliberately don't break and use the last valid
248 // reflection we find.
249 }
250 }
251 if (foundSatellite)
252 return std::make_tuple(indexedIntHKL, indexedMNP, indexedIntHKL.hklError());
253 else
254 return boost::none;
255}
256
268CombinedIndexingStats indexPeaks(const std::vector<IPeak *> &peaks, DblMatrix ub, const double mainTolerance,
269 const bool roundHKLs, const bool optimizeUB,
270 const Prop::SatelliteIndexingArgs &satelliteArgs) {
271 const auto nPeaks = peaks.size();
272 std::vector<V3D> qSample(nPeaks);
273 std::generate(std::begin(qSample), std::end(qSample),
274 [&peaks, i = 0u]() mutable { return peaks[i++]->getQSampleFrame(); });
275 if (optimizeUB) {
276 ub = optimizeUBMatrix(ub, qSample, mainTolerance);
277 }
278
279 CombinedIndexingStats stats;
280 ub.Invert();
281 for (auto i = 0u; i < peaks.size(); ++i) {
282 const auto peak = peaks[i];
283 V3D nominalHKL = IndexingUtils::CalculateMillerIndices(ub, qSample[i]);
284 if (IndexingUtils::ValidIndex(nominalHKL, mainTolerance)) {
285 stats.main.numIndexed++;
286 stats.main.error += nominalHKL.hklError() / 3.0;
287 if (roundHKLs) {
288 IndexingUtils::RoundHKL(nominalHKL);
289 }
290 peak->setHKL(nominalHKL);
291 peak->setIntHKL(nominalHKL);
292 peak->setIntMNP(V3D(0, 0, 0));
293 } else if (satelliteArgs.maxOrder > 0) {
294 auto result = indexSatellite(nominalHKL, satelliteArgs.maxOrder, satelliteArgs.modVectors,
295 satelliteArgs.tolerance, satelliteArgs.crossTerms);
296 if (result) {
297 const auto &satelliteInfo = result.get();
298 if (roundHKLs)
299 IndexingUtils::RoundHKL(nominalHKL);
300 peak->setHKL(nominalHKL);
301 peak->setIntHKL(std::get<0>(satelliteInfo));
302 peak->setIntMNP(std::get<1>(satelliteInfo));
303 stats.satellites.numIndexed++;
304 stats.satellites.error += std::get<2>(satelliteInfo) / 3.;
305 } else {
306 // clear these to make sure leftover values from previous index peaks
307 // run are not used
308 peak->setHKL(V3D(0, 0, 0));
309 peak->setIntHKL(V3D(0, 0, 0));
310 peak->setIntMNP(V3D(0, 0, 0));
311 }
312 } else {
313 peak->setHKL(V3D(0, 0, 0));
314 peak->setIntHKL(V3D(0, 0, 0));
315 peak->setIntMNP(V3D(0, 0, 0));
316 }
317 }
318 return stats;
319}
320
330void logIndexingResults(std::ostream &out, const CombinedIndexingStats &indexingInfo, const int runNo,
331 const size_t nPeaksTotal, const Prop::IndexPeaksArgs &args) {
332 if (runNo >= 0)
333 out << "Run " << runNo;
334 else
335 out << "All runs";
336 out << " indexed " << indexingInfo.totalNumIndexed() << " peaks out of " << nPeaksTotal;
337 if (args.satellites.maxOrder > 0) {
338 out << " of which, " << indexingInfo.main.numIndexed << " main Bragg peaks are indexed with tolerance of "
339 << args.mainTolerance << ", " << indexingInfo.satellites.numIndexed
340 << " satellite peaks are indexed with tolerance of " << args.satellites.tolerance << '\n';
341 out << " Average error in h,k,l for indexed peaks = " << indexingInfo.averageError() << '\n';
342 out << " Average error in h,k,l for indexed main peaks = " << indexingInfo.main.error << '\n';
343 out << " Average error in h,k,l for indexed satellite peaks = " << indexingInfo.satellites.error << '\n';
344 } else {
345 out << " with tolerance of " << args.mainTolerance << '\n';
346 out << " Average error in h,k,l for indexed peaks = " << indexingInfo.mainError() << '\n';
347 }
348}
349
350} // namespace
351
359 using Kernel::Direction;
360
361 // -- inputs --
362 this->declareProperty(std::make_unique<WorkspaceProperty<IPeaksWorkspace_sptr::element_type>>(Prop::PEAKSWORKSPACE,
363 "", Direction::InOut),
364 "Input Peaks Workspace");
365
366 auto mustBePositiveDbl = std::make_shared<BoundedValidator<double>>();
367 mustBePositiveDbl->setLower(0.0);
368 this->declareProperty(Prop::TOLERANCE, 0.15, mustBePositiveDbl, "Main peak indexing tolerance", Direction::Input);
369 this->declareProperty(Prop::SATE_TOLERANCE, 0.15, mustBePositiveDbl, "Satellite peak indexing tolerance",
371 this->declareProperty(Prop::ROUNDHKLS, true, "Round H, K and L values to integers");
372 this->declareProperty(Prop::COMMONUB, false, "Index all orientations with a common UB");
374 this->declareProperty(Prop::SAVEMODINFO, false,
375 "If true, update the OrientedLattice with the maxOrder, "
376 "modulation vectors & cross terms values input to the algorithm");
377
378 // -- outputs --
379 this->declareProperty(Prop::NUM_INDEXED, 0, "Gets set with the number of indexed peaks.", Direction::Output);
380 this->declareProperty(Prop::AVERAGE_ERR, 0.0, "Gets set with the average HKL indexing error.", Direction::Output);
381 this->declareProperty(Prop::MAIN_NUM_INDEXED, 0, "Gets set with the number of indexed main peaks.",
383 this->declareProperty(Prop::SATE_NUM_INDEXED, 0, "Gets set with the number of indexed main peaks.",
385 this->declareProperty(Prop::MAIN_ERR, 0.0, "Gets set with the average HKL indexing error of Main Peaks.",
387 this->declareProperty(Prop::SATE_ERR, 0.0, "Gets set with the average HKL indexing error of Satellite Peaks.",
389}
390
395std::map<std::string, std::string> IndexPeaks::validateInputs() {
396 std::map<std::string, std::string> helpMsgs;
397
398 IPeaksWorkspace_sptr ws = this->getProperty(Prop::PEAKSWORKSPACE);
399 try {
400 ws->sample().getOrientedLattice();
401 } catch (std::runtime_error &) {
402 helpMsgs[Prop::PEAKSWORKSPACE] = "No UB Matrix defined in the lattice.";
403 return helpMsgs;
404 }
405
406 const auto args = Prop::IndexPeaksArgs::parse(*this);
407 const bool isSave = this->getProperty(Prop::SAVEMODINFO);
408 const bool isMOZero = (args.satellites.maxOrder == 0);
409 bool isAllVecZero = true;
410 // parse() validates all the mod vectors. There should not be any modulated
411 // vector in modVectors is equal to (0, 0, 0)
412 for (size_t vecNo = 0; vecNo < args.satellites.modVectors.size(); vecNo++) {
413 if (args.satellites.modVectors[vecNo] != V3D(0.0, 0.0, 0.0)) {
414 isAllVecZero = false;
415 } else {
416 g_log.warning() << "Mod vector " << vecNo << " is invalid (0, 0, 0)"
417 << "\n";
418 }
419 }
420 if (isMOZero && !isAllVecZero) {
421 helpMsgs["MaxOrder"] = "Max Order cannot be zero if a Modulation Vector has been supplied.";
422 }
423 if (!isMOZero && isAllVecZero) {
424 helpMsgs["ModVector1"] = "At least one Modulation Vector must be supplied if Max Order set.";
425 }
426 if (isSave && isAllVecZero) {
427 helpMsgs[Prop::SAVEMODINFO] = "Modulation info cannot be saved with no "
428 "valid Modulation Vectors supplied.";
429 }
430 if (isSave && isMOZero) {
431 helpMsgs["MaxOrder"] = "Modulation info cannot be saved with Max Order = 0.";
432 }
433 return helpMsgs;
434}
435
440 const auto args = Prop::IndexPeaksArgs::parse(*this);
441
442 // quick exit
443 if (args.workspace->getNumberPeaks() == 0) {
444 g_log.warning("Empty peaks workspace. Nothing to index");
445 return;
446 }
447
448 // save modulation input if asked
449 if (args.storeModulationInfo) {
450 auto &lattice = args.workspace->mutableSample().getOrientedLattice();
451 lattice.setMaxOrder(args.satellites.maxOrder);
452 lattice.setCrossTerm(args.satellites.crossTerms);
453
454 if (args.satellites.modVectors.size() >= 1) {
455 lattice.setModVec1(args.satellites.modVectors[0]);
456 } else {
457 g_log.warning("empty modVector 1, skipping saving");
458 lattice.setModVec1(V3D(0.0, 0.0, 0.0));
459 }
460
461 if (args.satellites.modVectors.size() >= 2) {
462 lattice.setModVec2(args.satellites.modVectors[1]);
463 } else {
464 g_log.warning("empty modVector 2, skipping saving");
465 lattice.setModVec2(V3D(0.0, 0.0, 0.0));
466 }
467
468 if (args.satellites.modVectors.size() >= 3) {
469 lattice.setModVec3(args.satellites.modVectors[2]);
470 } else {
471 g_log.warning("empty modVector 3, skipping saving");
472 lattice.setModVec3(V3D(0.0, 0.0, 0.0));
473 }
474 }
475
476 CombinedIndexingStats indexingInfo;
477 const auto &lattice = args.workspace->sample().getOrientedLattice();
478 const auto &sampleUB = lattice.getUB();
479 if (args.commonUB) {
480 // Use sample UB an all peaks regardless of run
481 std::vector<IPeak *> allPeaksRef;
482 allPeaksRef.reserve(args.workspace->getNumberPeaks());
483 for (int i = 0; i < args.workspace->getNumberPeaks(); i++) {
484 allPeaksRef.emplace_back(args.workspace->getPeakPtr(i));
485 }
486 const bool optimizeUB{false};
487 indexingInfo = indexPeaks(allPeaksRef, sampleUB, args.mainTolerance, args.roundHKLs, optimizeUB, args.satellites);
488 } else {
489 // Use a UB optimized for each run
490 std::unordered_map<int, std::vector<IPeak *>> peaksPerRun;
491 for (int i = 0; i < args.workspace->getNumberPeaks(); i++)
492 peaksPerRun[args.workspace->getPeak(i).getRunNumber()].emplace_back(args.workspace->getPeakPtr(i));
493 if (peaksPerRun.size() < 2) {
494 g_log.warning("Peaks from only one run exist but CommonUBForAll=True so peaks will be indexed with an optimised "
495 "UB which will not be saved in the workspace.");
496 }
497 const bool optimizeUB{true};
498 for (const auto &runPeaks : peaksPerRun) {
499 const auto &peaks = runPeaks.second;
500 const auto indexedInRun =
501 indexPeaks(peaks, sampleUB, args.mainTolerance, args.roundHKLs, optimizeUB, args.satellites);
502 logIndexingResults(g_log.notice(), indexedInRun, runPeaks.first, peaks.size(), args);
503 indexingInfo += indexedInRun;
504 }
505 }
506
507 setProperty("NumIndexed", indexingInfo.totalNumIndexed());
508 setProperty("MainNumIndexed", indexingInfo.main.numIndexed);
509 setProperty("SateNumIndexed", indexingInfo.satellites.numIndexed);
510 setProperty("AverageError", indexingInfo.averageError());
511 setProperty("MainError", indexingInfo.mainError());
512 setProperty("SatelliteError", indexingInfo.satelliteError());
513
514 // Final results
515 logIndexingResults(g_log.notice(), indexingInfo, -1, args.workspace->getNumberPeaks(), args);
516 // Show the lattice parameters
517 g_log.notice() << args.workspace->sample().getOrientedLattice() << "\n";
518}
519
520} // namespace Mantid::Crystal
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
const std::vector< double > & rhs
const int64_t TOLERANCE(1000000)
const bool commonUB
Definition: IndexPeaks.cpp:117
double error
Definition: IndexPeaks.cpp:133
PeakIndexingStats main
Definition: IndexPeaks.cpp:164
int numIndexed
Definition: IndexPeaks.cpp:132
const std::vector< V3D > modVectors
Definition: IndexPeaks.cpp:56
const bool crossTerms
Definition: IndexPeaks.cpp:57
const bool storeModulationInfo
Definition: IndexPeaks.cpp:118
const bool roundHKLs
Definition: IndexPeaks.cpp:116
const int maxOrder
Definition: IndexPeaks.cpp:55
const SatelliteIndexingArgs satellites
Definition: IndexPeaks.cpp:119
IPeaksWorkspace_sptr workspace
Definition: IndexPeaks.cpp:114
const double mainTolerance
Definition: IndexPeaks.cpp:115
double tolerance
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
A property class for workspaces.
void init() override
Initialize the algorithm's properties.
Definition: IndexPeaks.cpp:354
std::map< std::string, std::string > validateInputs() override
Validate all inputs once set.
Definition: IndexPeaks.cpp:395
void exec() override
Execute the algorithm.
Definition: IndexPeaks.cpp:439
static void RoundHKL(Kernel::V3D &hkl)
Round all the components of a HKL objects to the nearest integer.
static void RoundHKLs(std::vector< Kernel::V3D > &hkl_list)
Round all the components of a list of V3D objects, to the nearest integer.
static double Optimize_UB(Kernel::DblMatrix &UB, const std::vector< Kernel::V3D > &hkl_vectors, const std::vector< Kernel::V3D > &q_vectors, std::vector< double > &sigabc)
Find the UB matrix that most nearly maps hkl to qxyz for 3 or more peaks.
static bool ValidIndex(const Kernel::V3D &hkl, double tolerance)
Check is hkl is within tolerance of integer (h,k,l) non-zero values.
static int CalculateMillerIndices(const Kernel::DblMatrix &UB, const std::vector< Kernel::V3D > &q_vectors, double tolerance, std::vector< Kernel::V3D > &miller_indices, double &ave_error)
Given a UB, get list of Miller indices for specifed Qs and tolerance.
ArrayLenghtValidator : Validate length of an array property.
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
BoundedValidator is a validator that requires the values to be between upper or lower bounds,...
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void notice(const std::string &msg)
Logs at notice level.
Definition: Logger.cpp:95
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
Class for 3D vectors.
Definition: V3D.h:34
std::shared_ptr< IPeaksWorkspace > IPeaksWorkspace_sptr
shared pointer to Mantid::API::IPeaksWorkspace
MatrixWorkspace_sptr MANTID_API_DLL operator+=(const MatrixWorkspace_sptr &lhs, const MatrixWorkspace_sptr &rhs)
Adds two workspaces.
std::vector< Kernel::V3D > addModulationVectors(const std::vector< double > &modVector1, const std::vector< double > &modVector2, const std::vector< double > &modVector3)
Create a list of valid modulation vectors from the input.
std::vector< MNPOffset > generateOffsetVectors(const std::vector< Kernel::V3D > &modVectors, const int maxOrder, const bool crossTerms)
Calculate a list of HKL offsets from the given modulation vectors.
std::vector< Kernel::V3D > validModulationVectors(const std::vector< double > &modVector1, const std::vector< double > &modVector2, const std::vector< double > &modVector3)
Create a list of valid modulation vectors from the input.
Mantid::Kernel::Matrix< double > DblMatrix
Definition: Matrix.h:206
String constants for algorithm's properties.
static void appendTo(API::IAlgorithm *alg)
Append the common set of properties that relate to modulation vectors to the given algorithm.
Describes the direction (within an algorithm) of a Property.
Definition: Property.h:50
@ InOut
Both an input & output workspace.
Definition: Property.h:55
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54