Mantid
Loading...
Searching...
No Matches
TransformHKL.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
19namespace Mantid::Crystal {
20// Register the algorithm into the AlgorithmFactory
21DECLARE_ALGORITHM(TransformHKL)
22
23using namespace Mantid::Kernel;
24using namespace Mantid::API;
25using namespace Mantid::DataObjects;
26using namespace Mantid::Geometry;
27
28const std::string TransformHKL::name() const { return "TransformHKL"; }
29
30int TransformHKL::version() const { return 1; }
31
32const std::string TransformHKL::category() const { return "Crystal\\Peaks"; }
33
37 this->declareProperty(std::make_unique<WorkspaceProperty<IPeaksWorkspace>>("PeaksWorkspace", "", Direction::InOut),
38 "Input Peaks Workspace");
39
40 auto mustBePositive = std::make_shared<BoundedValidator<double>>();
41 mustBePositive->setLower(0.0);
42
43 this->declareProperty(
44 std::make_unique<PropertyWithValue<double>>("Tolerance", 0.15, std::move(mustBePositive), Direction::Input),
45 "Indexing Tolerance (0.15)");
46
47 std::vector<double> identity_matrix(9, 0.0);
48 identity_matrix[0] = 1;
49 identity_matrix[4] = 1;
50 identity_matrix[8] = 1;
51 // clang-format off
52 auto threeBythree = std::make_shared<ArrayLengthValidator<double> >(9);
53 // clang-format on
54 this->declareProperty(
55 std::make_unique<ArrayProperty<double>>("HKLTransform", std::move(identity_matrix), std::move(threeBythree)),
56 "Specify 3x3 HKL transform matrix as a comma separated list of 9 "
57 "numbers");
58
59 this->declareProperty("FindError", true,
60 "Whether to obtain the error in the lattice parameters. "
61 "Set this to false if there are not enough peaks to do an optimization");
62
63 this->declareProperty(std::make_unique<PropertyWithValue<int>>("NumIndexed", 0, Direction::Output),
64 "Gets set with the number of indexed peaks.");
65
66 this->declareProperty(std::make_unique<PropertyWithValue<double>>("AverageError", 0.0, Direction::Output),
67 "Gets set with the average HKL indexing error.");
68}
69
73 IPeaksWorkspace_sptr ws = this->getProperty("PeaksWorkspace");
74 if (!ws) {
75 throw std::runtime_error("Could not read the peaks workspace");
76 }
77
78 OrientedLattice o_lattice = ws->mutableSample().getOrientedLattice();
79 Matrix<double> UB = o_lattice.getUB();
80
81 if (!IndexingUtils::CheckUB(UB)) {
82 throw std::runtime_error("ERROR: The stored UB is not a valid orientation matrix");
83 }
84
85 std::vector<double> tran_vec = getProperty("HKLTransform");
86 DblMatrix hkl_tran(tran_vec);
87
88 std::ostringstream str_stream;
89 str_stream << hkl_tran;
90 std::string hkl_tran_string = str_stream.str();
91 g_log.notice() << "Applying Tranformation " << hkl_tran_string << '\n';
92
93 if (hkl_tran.numRows() != 3 || hkl_tran.numCols() != 3) {
94 throw std::runtime_error("ERROR: The specified transform must be a 3 X 3 matrix.\n" + hkl_tran_string);
95 }
96
97 DblMatrix hkl_tran_inverse(hkl_tran);
98 double det = hkl_tran_inverse.Invert();
99
100 if (fabs(det) < 1.0e-5) {
101 throw std::runtime_error("ERROR: The specified matrix is invalid (essentially singular.)" + hkl_tran_string);
102 }
103 if (det < 0) {
104 std::ostringstream error_stream;
105 error_stream << hkl_tran;
106 throw std::runtime_error("ERROR: The determinant of the matrix is negative.\n" + str_stream.str());
107 }
108 double tolerance = this->getProperty("Tolerance");
109
110 // Transform looks OK so update UB and
111 // transform the hkls
112 UB = UB * hkl_tran_inverse;
113 g_log.notice() << "Transformed UB = " << UB << '\n';
114 o_lattice.setUB(UB);
115
116 Matrix<double> modHKL = o_lattice.getModHKL();
117
118 Matrix<double> modUB = UB * modHKL;
119
120 o_lattice.setModHKL(hkl_tran * modHKL);
121
122 std::vector<double> sigabc(6);
123
124 bool redetermine_error = this->getProperty("FindError");
125
126 if (redetermine_error) {
127 SelectCellWithForm::DetermineErrors(sigabc, UB, modUB, ws, tolerance);
128 o_lattice.setError(sigabc[0], sigabc[1], sigabc[2], sigabc[3], sigabc[4], sigabc[5]);
129 } else {
130 o_lattice.setError(0, 0, 0, 0, 0, 0);
131 }
132
133 ws->mutableSample().setOrientedLattice(std::make_unique<OrientedLattice>(o_lattice));
134
135 int n_peaks = ws->getNumberPeaks();
136
137 // transform all the HKLs and record the new HKL
138 // and q-vectors for peaks ORIGINALLY indexed
139 int num_indexed = 0;
140 std::vector<V3D> miller_indices;
141 std::vector<V3D> q_vectors;
142 for (int i = 0; i < n_peaks; i++) {
143 IPeak &peak = ws->getPeak(i);
144 V3D hkl(peak.getHKL());
145 V3D ihkl(peak.getIntHKL());
146 peak.setIntHKL(hkl_tran * ihkl);
147 miller_indices.emplace_back(hkl_tran * hkl - modHKL * peak.getIntMNP());
148 peak.setHKL(hkl_tran * hkl);
149 q_vectors.emplace_back(peak.getQSampleFrame() - modUB * peak.getIntMNP() * 2 * M_PI);
150 num_indexed++;
151 }
152
153 double average_error = IndexingUtils::IndexingError(UB, miller_indices, q_vectors);
154
155 // Tell the user what happened.
156 g_log.notice() << o_lattice << "\n";
157 g_log.notice() << "Transformed Miller indices on previously valid indexed Peaks.\n";
158 g_log.notice() << "Set hkl to 0,0,0 on peaks previously indexed out of tolerance.\n";
159 g_log.notice() << "Now, " << num_indexed << " are indexed with average error " << average_error << "\n";
160
161 // Save output properties
162 this->setProperty("NumIndexed", num_indexed);
163 this->setProperty("AverageError", average_error);
164}
165
166} // namespace Mantid::Crystal
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
#define fabs(x)
Definition: Matrix.cpp:22
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.
static Kernel::Matrix< double > DetermineErrors(std::vector< double > &sigabc, const Kernel::Matrix< double > &UB, const Kernel::Matrix< double > &modUB, const API::IPeaksWorkspace_sptr &ws, double tolerance)
const std::string name() const override
Algorithm's name for identification.
void init() override
Initialise the properties.
void exec() override
Run the algorithm.
const std::string category() const override
Algorithm's category for identification.
int version() const override
Algorithm's version for identification.
Structure describing a single-crystal peak.
Definition: IPeak.h:26
virtual void setHKL(double H, double K, double L)=0
virtual Mantid::Kernel::V3D getIntHKL() const =0
virtual Mantid::Kernel::V3D getQSampleFrame() const =0
virtual void setIntHKL(const Mantid::Kernel::V3D &HKL)=0
virtual Mantid::Kernel::V3D getHKL() const =0
virtual Mantid::Kernel::V3D getIntMNP() const =0
static double IndexingError(const Kernel::DblMatrix &UB, const std::vector< Kernel::V3D > &hkls, const std::vector< Kernel::V3D > &q_vectors)
Find the average indexing error for UB with the specified q's and hkls.
static bool CheckUB(const Kernel::DblMatrix &UB)
Check that the specified UB is reasonable for an orientation matrix.
Class to implement UB matrix.
void setUB(const Kernel::DblMatrix &newUB)
Sets the UB matrix and recalculates lattice parameters.
const Kernel::DblMatrix & getUB() const
Get the UB matrix.
const Kernel::DblMatrix & getModHKL() const
Get modulation vectors for satellites.
Definition: UnitCell.cpp:548
void setModHKL(double _dh1, double _dk1, double _dl1, double _dh2, double _dk2, double _dl2, double _dh3, double _dk3, double _dl3)
Set modulation vectors for satellites.
Definition: UnitCell.cpp:353
void setError(double _aerr, double _berr, double _cerr, double _alphaerr, double _betaerr, double _gammaerr, const int angleunit=angDegrees)
Set lattice parameter errors.
Definition: UnitCell.cpp:325
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
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
Numerical Matrix class.
Definition: Matrix.h:42
T Invert()
LU inversion routine.
Definition: Matrix.cpp:924
size_t numRows() const
Return the number of rows in the matrix.
Definition: Matrix.h:144
size_t numCols() const
Return the number of columns in the matrix.
Definition: Matrix.h:147
std::string str() const
Convert the matrix into a simple linear string expression.
Definition: Matrix.cpp:1564
The concrete, templated class for properties.
Class for 3D vectors.
Definition: V3D.h:34
std::shared_ptr< IPeaksWorkspace > IPeaksWorkspace_sptr
shared pointer to Mantid::API::IPeaksWorkspace
@ 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