Mantid
Loading...
Searching...
No Matches
GoniometerAnglesFromPhiRotation.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 +
10#include "MantidAPI/IFunction.h"
11#include "MantidAPI/Sample.h"
17
19
20namespace Mantid::Crystal {
21
22// Register the algorithm into the AlgorithmFactory
23DECLARE_ALGORITHM(GoniometerAnglesFromPhiRotation)
24
25using namespace Mantid::Kernel;
26using namespace Mantid::API;
27using namespace Mantid::DataObjects;
28using namespace Mantid::Geometry;
29
31 declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace>>("PeaksWorkspace1", "", Kernel::Direction::Input),
32 "Input Peaks Workspace for Run 1");
33
34 declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace>>("PeaksWorkspace2", "", Kernel::Direction::InOut),
35 "Input Peaks Workspace for Run 2");
36
37 declareProperty("Tolerance", .12, "Integer offset for h,k,and l values to be considered valid.(def=.12)");
38
39 declareProperty("MIND", -1.0, "Minimium d-spacing to consider,(def=-1)");
40 declareProperty("MAXD", -1.0, "Maximum d-spacing to consider,(def=-1)");
41
42 declareProperty("Run1Phi", 0.0, "Phi for Run 1(def=0.0)");
43
44 declareProperty(std::string("Phi2"), 0.0, std::string("Phi angle for Run2(def=0.0)"), Kernel::Direction::InOut);
45
46 declareProperty("Chi2", 0.0, "Chi angle for Run2", Kernel::Direction::Output);
47 declareProperty("Omega2", 0.0, "Omega angle for Run2", Kernel::Direction::Output);
48
49 declareProperty("NIndexed", 0, "Number peaks indexed", Kernel::Direction::Output);
50
51 declareProperty("AvErrIndex", 0.0, "Average abs offset from integer values for indexed peaks",
53
54 declareProperty("AvErrAll", 0.0, "Average abs offset from integer values for all peaks", Kernel::Direction::Output);
55}
56
74 int &Nindexed, double &AvErrIndexed, double &AvErrorAll,
75 double tolerance) const {
76 Kernel::Matrix<double> InvUB2(UBraw);
77
78 InvUB2.Invert();
79 InvUB2 /= (2 * M_PI);
80
81 Nindexed = 0;
82 double TotOffsetIndx = 0;
83 double TotOffsetAll = 0;
84 int Npeaks = Peaks->getNumberPeaks();
85
86 for (int i = 0; i < Npeaks; i++) {
87 V3D hkl = InvUB2 * Peaks->getPeak(i).getQLabFrame();
88 double maxOffset = 0;
89 for (int k = 0; k < 3; k++) {
90 double offset = hkl[k] - floor(hkl[k]);
91 if (offset > .5)
92 offset -= 1;
93 offset = fabs(offset);
94 if (offset > maxOffset)
95 maxOffset = offset;
96 }
97
98 if (maxOffset < tolerance) {
99 Nindexed++;
100 TotOffsetIndx += maxOffset;
101 }
102
103 TotOffsetAll += maxOffset;
104 }
105
106 if (Nindexed > 0)
107
108 AvErrIndexed = TotOffsetIndx / Nindexed;
109
110 else
111
112 AvErrIndexed = -1.0;
113
114 if (Npeaks > 0)
115
116 AvErrorAll = TotOffsetAll / Npeaks;
117
118 else
119
120 AvErrorAll = -1.0;
121}
123
124 PeaksWorkspace_sptr PeaksRun1 = getProperty("PeaksWorkspace1");
125 PeaksWorkspace_sptr PeaksRun2 = getProperty("PeaksWorkspace2");
126
127 double Tolerance = getProperty("Tolerance");
128
129 Kernel::Matrix<double> Gon1(3, 3);
130 Kernel::Matrix<double> Gon2(3, 3);
131 if (!CheckForOneRun(PeaksRun1, Gon1) || !CheckForOneRun(PeaksRun2, Gon2)) {
132 g_log.error("Each peaks workspace MUST have only one run");
133 throw std::invalid_argument("Each peaks workspace MUST have only one run");
134 }
135
137
138 bool Run1HasOrientedLattice = true;
139 if (!PeaksRun1->sample().hasOrientedLattice()) {
140
141 Run1HasOrientedLattice = false;
142
143 const std::string fft("FindUBUsingFFT");
144 auto findUB = createChildAlgorithm(fft);
145 findUB->initialize();
146 findUB->setProperty<PeaksWorkspace_sptr>("PeaksWorkspace", getProperty("PeaksWorkspace1"));
147 findUB->setProperty("MIND", static_cast<double>(getProperty("MIND")));
148 findUB->setProperty("MAXD", static_cast<double>(getProperty("MAXD")));
149 findUB->setProperty("Tolerance", Tolerance);
150
151 findUB->executeAsChildAlg();
152
153 if (!PeaksRun1->sample().hasOrientedLattice()) {
154 g_log.notice(std::string("Could not find UB for ") + std::string(PeaksRun1->getName()));
155 throw std::invalid_argument(std::string("Could not find UB for ") + std::string(PeaksRun1->getName()));
156 }
157 }
158 //-------------get UB raw :No goniometer----------------
159
160 UB1 = PeaksRun1->sample().getOrientedLattice().getUB();
161
162 UB1 = getUBRaw(UB1, Gon1);
163
164 int N1;
165 double avErrIndx, avErrAll;
166 IndexRaw(PeaksRun1, UB1, N1, avErrIndx, avErrAll, Tolerance);
167
168 if (N1 < .6 * PeaksRun1->getNumberPeaks()) {
169 g_log.notice(std::string("UB did not index well for ") + std::string(PeaksRun1->getName()));
170 throw std::invalid_argument(std::string("UB did not index well for ") + std::string(PeaksRun1->getName()));
171 }
172
173 //----------------------------------------------
174
175 auto lat2 = std::make_unique<OrientedLattice>(PeaksRun1->sample().getOrientedLattice());
176 lat2->setUB(UB1);
177 PeaksRun2->mutableSample().setOrientedLattice(std::move(lat2));
178
179 if (!Run1HasOrientedLattice)
180 PeaksRun1->mutableSample().setOrientedLattice(nullptr);
181
182 double dphi = static_cast<double>(getProperty("Phi2")) - static_cast<double>(getProperty("Run1Phi"));
183 Kernel::Matrix<double> Gon22(3, 3, true);
184
185 for (int i = 0; i < PeaksRun2->getNumberPeaks(); i++) {
186 PeaksRun2->getPeak(i).setGoniometerMatrix(Gon22);
187 }
188
189 int RunNum = PeaksRun2->getPeak(0).getRunNumber();
190 std::string RunNumStr = std::to_string(RunNum);
191 int Npeaks = PeaksRun2->getNumberPeaks();
192
193 // n indexed, av err, phi, chi,omega
194 std::array<double, 5> MinData = {{0., 0., 0., 0., 0.}};
195 MinData[0] = 0.0;
196 std::vector<V3D> directionList = IndexingUtils::MakeHemisphereDirections(50);
197
199
200 for (auto dir : directionList)
201 for (int sgn = 1; sgn > -2; sgn -= 2) {
202 dir.normalize();
203 Quat Q(sgn * dphi, dir);
204 Q.normalize();
206
207 int Nindexed;
208 double dummyAvErrIndx, dummyAvErrAll;
209 IndexRaw(PeaksRun2, Rot * UB1, Nindexed, dummyAvErrIndx, dummyAvErrAll, Tolerance);
210
211 if (Nindexed > MinData[0]) {
212 MinData[0] = Nindexed;
213 MinData[1] = sgn;
214 MinData[2] = dir[0];
215
216 MinData[3] = dir[1];
217 MinData[4] = dir[2];
218 }
219 }
220
221 g_log.debug() << "Best direction unOptimized is (" << (MinData[1] * MinData[2]) << "," << (MinData[1] * MinData[3])
222 << "," << (MinData[1] * MinData[4]) << ")\n";
223
224 //----------------------- Optimize around best----------------------------
225
226 auto ws = createWorkspace<Workspace2D>(1, 3 * Npeaks, 3 * Npeaks);
227
228 MantidVec Xvals;
229
230 for (int i = 0; i < Npeaks; ++i) {
231 Xvals.emplace_back(i);
232 Xvals.emplace_back(i);
233 Xvals.emplace_back(i);
234 }
235
236 ws->setPoints(0, Xvals);
237
238 //--------------------Set up other Fit function arguments------------------
239 V3D dir(MinData[2], MinData[3], MinData[4]);
240 dir.normalize();
241 Quat Q(MinData[1] * dphi, dir);
242 Q.normalize();
244
245 Goniometer Gon(Rot);
246 std::vector<double> omchiphi = Gon.getEulerAngles("yzy");
247 MinData[2] = omchiphi[2];
248 MinData[3] = omchiphi[1];
249 MinData[4] = omchiphi[0];
250
251 std::string FunctionArgs = "name=PeakHKLErrors, PeakWorkspaceName=" + PeaksRun2->getName() + ",OptRuns=" + RunNumStr +
252 ",phi" + RunNumStr + "=" + boost::lexical_cast<std::string>(MinData[2]) + ",chi" +
253 RunNumStr + "=" + boost::lexical_cast<std::string>(MinData[3]) + ",omega" + RunNumStr +
254 "=" + boost::lexical_cast<std::string>(MinData[4]);
255
256 std::string Constr = boost::lexical_cast<std::string>(MinData[2] - 5) + "<phi" + RunNumStr + "<" +
257 boost::lexical_cast<std::string>(MinData[2] + 5);
258 Constr += "," + boost::lexical_cast<std::string>(MinData[3] - 5) + "<chi" + RunNumStr + "<" +
259 boost::lexical_cast<std::string>(MinData[3] + 5) + ",";
260
261 Constr += boost::lexical_cast<std::string>(MinData[4] - 5) + "<omega" + RunNumStr + "<" +
262 boost::lexical_cast<std::string>(MinData[4] + 5);
263
264 std::string Ties = "SampleXOffset=0.0,SampleYOffset=0.0,SampleZOffset=0.0,"
265 "GonRotx=0.0,GonRoty=0.0,GonRotz=0.0";
266
267 std::shared_ptr<Algorithm> Fit = createChildAlgorithm("Fit");
268
269 Fit->initialize();
270 Fit->setProperty("Function", FunctionArgs);
271 Fit->setProperty("Ties", Ties);
272 Fit->setProperty("Constraints", Constr);
273 Fit->setProperty("InputWorkspace", ws);
274 Fit->setProperty("CreateOutput", true);
275
276 std::string outputName = "out";
277
278 Fit->setProperty("Output", outputName);
279
281
282 std::shared_ptr<API::ITableWorkspace> results = Fit->getProperty("OutputParameters");
283 double chisq = Fit->getProperty("OutputChi2overDoF");
284
285 MinData[0] = chisq;
286 MinData[2] = results->Double(6, 1);
287 MinData[3] = results->Double(7, 1);
288 MinData[4] = results->Double(8, 1);
289
290 g_log.debug() << "Best direction Optimized is (" << (MinData[2]) << "," << (MinData[3]) << "," << (MinData[4])
291 << ")\n";
292
293 // ---------------------Find number indexed -----------------------
294 Quat Q1 = Quat(MinData[4], V3D(0, 1, 0)) * Quat(MinData[3], V3D(0, 0, 1)) * Quat(MinData[2], V3D(0, 1, 0));
295
296 int Nindexed;
298 IndexRaw(PeaksRun2, Mk * UB1, Nindexed, avErrIndx, avErrAll, Tolerance);
299
300 //------------------------------------ Convert/Save Results
301 //-----------------------------
302
303 double deg, ax1, ax2, ax3;
304 Q1.getAngleAxis(deg, ax1, ax2, ax3);
305 if (dphi * deg < 0) {
306 deg = -deg;
307 ax1 = -ax1;
308 ax2 = -ax2;
309 ax3 = -ax3;
310 }
311
312 double phi2 = static_cast<double>(getProperty("Run1Phi")) + dphi;
313 double chi2 = acos(ax2) / M_PI * 180;
314 double omega2 = atan2(ax3, -ax1) / M_PI * 180;
315
316 g_log.notice() << "============================ Results ============================\n";
317 g_log.notice() << " phi,chi, and omega= (" << phi2 << "," << chi2 << "," << omega2 << ")\n";
318 g_log.notice() << " #indexed =" << Nindexed << '\n';
319 g_log.notice() << " ==============================================\n";
320
321 setProperty("Phi2", phi2);
322 setProperty("Chi2", chi2);
323 setProperty("Omega2", omega2);
324
325 setProperty("NIndexed", Nindexed);
326
327 setProperty("AvErrIndex", avErrIndx);
328
329 setProperty("AvErrAll", avErrAll);
330
331 Q1 = Quat(omega2, V3D(0, 1, 0)) * Quat(chi2, V3D(0, 0, 1)) * Quat(phi2, V3D(0, 1, 0));
333 for (int i = 0; i < PeaksRun2->getNumberPeaks(); i++) {
334 PeaksRun2->getPeak(i).setGoniometerMatrix(Gon2a);
335 }
336
337 auto latt2 = std::make_unique<OrientedLattice>(PeaksRun2->mutableSample().getOrientedLattice());
338 Rot.Invert();
339 Gon2a.Invert();
340 latt2->setUB(Gon2a * Mk * UB1);
341
342 PeaksRun2->mutableSample().setOrientedLattice(std::move(latt2));
343}
344
352 Kernel::Matrix<double> &GoniometerMatrix) const {
353
354 int RunNumber = -1;
355 for (int peak = 0; peak < Peaks->getNumberPeaks(); peak++) {
356 int thisRunNum = Peaks->getPeak(peak).getRunNumber();
357 GoniometerMatrix = Peaks->getPeak(peak).getGoniometerMatrix();
358
359 if (RunNumber < 0)
360
361 RunNumber = thisRunNum;
362
363 else if (thisRunNum != RunNumber)
364
365 return false;
366 }
367
368 return true;
369}
370
379 const Kernel::Matrix<double> &GoniometerMatrix) const {
380 return GoniometerMatrix * UB;
381}
382
383} // 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
void initialize() override
Initialization method invoked by the framework.
Definition: Algorithm.cpp:283
virtual std::shared_ptr< Algorithm > createChildAlgorithm(const std::string &name, const double startProgress=-1., const double endProgress=-1., const bool enableLogging=true, const int &version=-1)
Create a Child Algorithm.
Definition: Algorithm.cpp:842
Kernel::Logger & g_log
Definition: Algorithm.h:451
void executeAsChildAlg() override
Execute as a Child Algorithm.
Definition: Algorithm.cpp:769
A property class for workspaces.
bool CheckForOneRun(const DataObjects::PeaksWorkspace_sptr &Peaks, Kernel::Matrix< double > &GoniometerMatrix) const
Checks that a PeaksWorkspace has only one run.
void IndexRaw(const DataObjects::PeaksWorkspace_sptr &Peaks, const Kernel::Matrix< double > &UBraw, int &Nindexed, double &AvErrIndexed, double &AvErrorAll, double tolerance) const
Calculate indexing stats if the peaks had been indexed with given UBraw by NOT applying the goniomete...
Kernel::Matrix< double > getUBRaw(const Kernel::Matrix< double > &UB, const Kernel::Matrix< double > &GoniometerMatrix) const
Returns the raw UB, its inverse indexes using Q lab.
A generic fitting algorithm.
Definition: Fit.h:78
Class to represent a particular goniometer setting, which is described by the rotation matrix.
Definition: Goniometer.h:55
std::vector< double > getEulerAngles(const std::string &convention="YZX")
Return Euler angles acording to a convention.
Definition: Goniometer.cpp:282
static std::vector< Kernel::V3D > MakeHemisphereDirections(int n_steps)
Make list of direction vectors uniformly distributed over a hemisphere.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
void notice(const std::string &msg)
Logs at notice level.
Definition: Logger.cpp:95
void error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
Numerical Matrix class.
Definition: Matrix.h:42
T Invert()
LU inversion routine.
Definition: Matrix.cpp:924
Class for quaternions.
Definition: Quat.h:39
void normalize()
Normalize.
Definition: Quat.cpp:340
void getAngleAxis(double &_deg, double &_ax0, double &_ax1, double &ax2) const
Extracts the angle of roatation and axis.
Definition: Quat.cpp:135
std::vector< double > getRotation(bool check_normalisation=false, bool throw_on_errors=false) const
returns the rotation matrix defined by this quaternion as an 9-point
Definition: Quat.cpp:453
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
Class for 3D vectors.
Definition: V3D.h:34
double normalize()
Make a normalized vector (return norm value)
Definition: V3D.cpp:130
std::shared_ptr< PeaksWorkspace > PeaksWorkspace_sptr
Typedef for a shared pointer to a peaks workspace.
constexpr double Tolerance
Standard tolerance value.
Definition: Tolerance.h:12
std::vector< double > MantidVec
typedef for the data storage used in Mantid matrix workspaces
Definition: cow_ptr.h:172
std::string to_string(const wide_integer< Bits, Signed > &n)
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