Mantid
Loading...
Searching...
No Matches
FindUBUsingFFT.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"
15
16namespace Mantid::Crystal {
17// Register the algorithm into the AlgorithmFactory
18DECLARE_ALGORITHM(FindUBUsingFFT)
19
20using namespace Mantid::Kernel;
21using namespace Mantid::API;
22using namespace Mantid::Geometry;
23
24const std::string FindUBUsingFFT::name() const { return "FindUBUsingFFT"; }
25
26int FindUBUsingFFT::version() const { return 1; }
27
28const std::string FindUBUsingFFT::category() const { return "Crystal\\UBMatrix"; }
29
33 this->declareProperty(std::make_unique<WorkspaceProperty<IPeaksWorkspace>>("PeaksWorkspace", "", Direction::InOut),
34 "Input Peaks Workspace");
35
36 auto mustBePositive = std::make_shared<BoundedValidator<double>>();
37 mustBePositive->setLower(0.0);
38
39 // use negative values, force user to input all parameters
40 this->declareProperty("MinD", -1.0, mustBePositive, "Lower Bound on Lattice Parameters a, b, c");
41 this->declareProperty("MaxD", -1.0, mustBePositive, "Upper Bound on Lattice Parameters a, b, c");
42 this->declareProperty("Tolerance", 0.15, mustBePositive, "Indexing Tolerance (0.15)");
43 this->declareProperty("Iterations", 4, "Iterations to refine UB (4)");
44 this->declareProperty("DegreesPerStep", 1.5,
45 "The resolution of the search through possible "
46 "orientations is specified by this parameter. One to "
47 "two degrees per step is usually adequate.");
48}
49
53 double min_d = this->getProperty("MinD");
54 double max_d = this->getProperty("MaxD");
55 double tolerance = this->getProperty("Tolerance");
56 int iterations = this->getProperty("Iterations");
57
58 double degrees_per_step = this->getProperty("DegreesPerStep");
59
60 IPeaksWorkspace_sptr ws = this->getProperty("PeaksWorkspace");
61
62 int n_peaks = ws->getNumberPeaks();
63
64 std::vector<V3D> q_vectors;
65 for (int i = 0; i < n_peaks; i++) {
66 q_vectors.emplace_back(ws->getPeak(i).getQSampleFrame());
67 }
68
69 Matrix<double> UB(3, 3, false);
70 double error = IndexingUtils::Find_UB(UB, q_vectors, min_d, max_d, tolerance, degrees_per_step, iterations);
71
72 g_log.notice() << "Error = " << error << '\n';
73 g_log.notice() << "UB = " << UB << '\n';
74
75 if (!IndexingUtils::CheckUB(UB)) // UB not found correctly
76 {
77 g_log.notice(std::string("Found Invalid UB...peaks used might not be linearly independent"));
78 g_log.notice(std::string("UB NOT SAVED."));
79 } else // tell user how many would be indexed
80 {
81 std::vector<double> sigabc(7);
82 std::vector<V3D> miller_ind;
83 std::vector<V3D> indexed_qs;
84 double fit_error;
85 miller_ind.reserve(q_vectors.size());
86 indexed_qs.reserve(q_vectors.size());
87 IndexingUtils::GetIndexedPeaks(UB, q_vectors, tolerance, miller_ind, indexed_qs, fit_error);
88 IndexingUtils::Optimize_UB(UB, miller_ind, indexed_qs, sigabc);
89
90 int num_indexed = IndexingUtils::NumberIndexed(UB, q_vectors, tolerance);
91 g_log.notice() << "New UB will index " << num_indexed << " Peaks out of " << n_peaks << " with tolerance of "
92 << std::setprecision(3) << std::setw(5) << tolerance << "\n";
93
94 auto lattice = std::make_unique<OrientedLattice>();
95 lattice->setUB(UB);
96 lattice->setError(sigabc[0], sigabc[1], sigabc[2], sigabc[3], sigabc[4], sigabc[5]);
97 // Show the modified lattice parameters
98 g_log.notice() << *lattice << "\n";
99 ws->mutableSample().setOrientedLattice(std::move(lattice));
100 }
101}
102
103} // namespace Mantid::Crystal
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double error
Definition: IndexPeaks.cpp:133
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.
int version() const override
Algorithm's version for identification.
void exec() override
Run the algorithm.
const std::string category() const override
Algorithm's category for identification.
const std::string name() const override
Algorithm's name for identification.
void init() override
Initialise the properties.
static int GetIndexedPeaks(const Kernel::DblMatrix &UB, const std::vector< Kernel::V3D > &q_vectors, double required_tolerance, std::vector< Kernel::V3D > &miller_indices, std::vector< Kernel::V3D > &indexed_qs, double &fit_error)
Get lists of indices and Qs for peaks indexed by the specified UB matrix.
static double Find_UB(Kernel::DblMatrix &UB, const std::vector< Kernel::V3D > &q_vectors, OrientedLattice &lattice, double required_tolerance, int base_index, size_t num_initial, double degrees_per_step, bool fixAll=false, int iterations=1)
Find the UB matrix that most nearly indexes the specified qxyz values given the lattice parameters.
static bool CheckUB(const Kernel::DblMatrix &UB)
Check that the specified UB is reasonable for an orientation matrix.
static int NumberIndexed(const Kernel::DblMatrix &UB, const std::vector< Kernel::V3D > &q_vectors, double tolerance)
Calculate the number of Q vectors that are mapped to integer indices by UB.
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.
void notice(const std::string &msg)
Logs at notice level.
Definition: Logger.cpp:95
Numerical Matrix class.
Definition: Matrix.h:42
std::shared_ptr< IPeaksWorkspace > IPeaksWorkspace_sptr
shared pointer to Mantid::API::IPeaksWorkspace
@ InOut
Both an input & output workspace.
Definition: Property.h:55