Mantid
Loading...
Searching...
No Matches
SaveIsawUB.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/Sample.h"
11
12#include <fstream>
13#include <iomanip>
14
17
18namespace Mantid::Crystal {
19
20// Register the algorithm into the AlgorithmFactory
21DECLARE_ALGORITHM(SaveIsawUB)
22
23using namespace Mantid::Kernel;
24using namespace Mantid::API;
25using namespace std;
26
30 declareProperty(std::make_unique<WorkspaceProperty<Workspace>>("InputWorkspace", "", Direction::Input),
31 "An input workspace containing the orientation matrix.");
32
33 const std::vector<std::string> exts{".mat", ".ub", ".txt"};
34 declareProperty(std::make_unique<FileProperty>("Filename", "", FileProperty::Save, exts),
35 "Path to an ISAW-style UB matrix text file.");
36}
37
39 double Volume;
40 double latticeParams[6] = {lattice.a(), lattice.b(), lattice.c(), lattice.alpha(), lattice.beta(), lattice.gamma()};
41 double lattice_errors[6] = {lattice.errora(), lattice.errorb(), lattice.errorc(),
42 lattice.erroralpha(), lattice.errorbeta(), lattice.errorgamma()};
43 if (lattice.volume() <= 0) {
44
45 double xA = cos(lattice.alpha() / 180. * M_PI);
46 double xB = cos(lattice.beta() / 180. * M_PI);
47 double xC = cos(lattice.gamma() / 180. * M_PI);
48 Volume = lattice.a() * lattice.b() * lattice.c() * sqrt(1 - xA * xA - xB * xB - xC * xC + 2 * xA * xB * xC);
49 } else
50 Volume = lattice.volume();
51
52 double dV = 0;
53 for (int i = 0; i < 3; i++) {
54 double U = (Volume / latticeParams[i] * lattice_errors[i]);
55 dV += U * U;
56 }
57
58 double U = (lattice_errors[3]) * (sin(2 * latticeParams[3] / 180. * M_PI) - sin(latticeParams[3] / 180. * M_PI) *
59 cos(latticeParams[4] / 180 * M_PI) *
60 cos(latticeParams[5] / 180 * M_PI));
61 dV += U * U;
62 U = (lattice_errors[4]) *
63 (sin(2 * latticeParams[4] / 180. * M_PI) -
64 sin(latticeParams[4] / 180. * M_PI) * cos(latticeParams[3] / 180 * M_PI) * cos(latticeParams[5] / 180 * M_PI));
65 dV += U * U;
66 U = (lattice_errors[5]) *
67 (sin(2 * latticeParams[5] / 180. * M_PI) -
68 sin(latticeParams[5] / 180. * M_PI) * cos(latticeParams[4] / 180 * M_PI) * cos(latticeParams[3] / 180 * M_PI));
69 dV += U * U;
70 dV = sqrt(dV);
71
72 return dV;
73}
74
78 try {
79 Workspace_sptr ws1 = getProperty("InputWorkspace");
81 MultipleExperimentInfos_sptr MDWS = std::dynamic_pointer_cast<MultipleExperimentInfos>(ws1);
82 if (MDWS != nullptr) {
83 ws = MDWS->getExperimentInfo(0);
84 } else {
85 ws = std::dynamic_pointer_cast<ExperimentInfo>(ws1);
86 }
87
88 if (!ws)
89 throw std::invalid_argument("Must specify either a MatrixWorkspace or a "
90 "PeaksWorkspace or a MDWorkspace.");
91
92 if (!ws->sample().hasOrientedLattice())
93 throw std::invalid_argument("Workspace must have an oriented lattice to save");
94
95 std::string Filename = getProperty("Filename");
96
97 ofstream out;
98 out.open(Filename.c_str());
99
100 OrientedLattice lattice = ws->sample().getOrientedLattice();
101 Kernel::DblMatrix ub = lattice.getUB();
102 Kernel::DblMatrix modub = lattice.getModUB();
103
104 // Write the ISAW UB matrix
105 const int beam = 2;
106 const int up = 1;
107 const int back = 0;
108 out << fixed;
109
110 for (size_t basis = 0; basis < 3; basis++) {
111 out << setw(11) << setprecision(8) << ub[beam][basis] << setw(12) << setprecision(8) << ub[back][basis]
112 << setw(12) << setprecision(8) << ub[up][basis] << " \n";
113 }
114
115 int ModDim = 0;
116 for (int i = 0; i < 3; i++) {
117 if (lattice.getModVec(i) == V3D(0, 0, 0))
118 continue;
119 else
120 ModDim++;
121 }
122
123 if (ModDim > 0) {
124 out << "ModUB: \n";
125 for (size_t basis = 0; basis < 3; basis++) {
126 out << setw(11) << setprecision(8) << modub[beam][basis] << setw(12) << setprecision(8) << modub[back][basis]
127 << setw(12) << setprecision(8) << modub[up][basis] << " \n";
128 }
129 }
130
131 // out << "Lattice Parameters: \n";
132 out << setw(11) << setprecision(4) << lattice.a() << setw(12) << setprecision(4) << lattice.b() << setw(12)
133 << setprecision(4) << lattice.c() << setw(12) << setprecision(4) << lattice.alpha() << setw(12)
134 << setprecision(4) << lattice.beta() << setw(12) << setprecision(4) << lattice.gamma() << setw(12)
135 << setprecision(4) << lattice.volume() << " \n";
136 double ErrorVolume = getErrorVolume(lattice);
137 out << setw(11) << setprecision(4) << lattice.errora() << setw(12) << setprecision(4) << lattice.errorb()
138 << setw(12) << setprecision(4) << lattice.errorc() << setw(12) << setprecision(4) << lattice.erroralpha()
139 << setw(12) << setprecision(4) << lattice.errorbeta() << setw(12) << setprecision(4) << lattice.errorgamma()
140 << setw(12) << setprecision(4) << ErrorVolume << " \n";
141
142 out << "\n";
143 if (ModDim >= 1) {
144 out << "Modulation Vector 1: " << setw(12) << setprecision(4) << lattice.getdh(0) << setw(12) << setprecision(4)
145 << lattice.getdk(0) << setw(12) << setprecision(4) << lattice.getdl(0) << " \n";
146
147 out << "Modulation Vector 1 error: " << setw(6) << setprecision(4) << lattice.getdherr(0) << setw(12)
148 << setprecision(4) << lattice.getdkerr(0) << setw(12) << setprecision(4) << lattice.getdlerr(0) << " \n";
149 }
150 if (ModDim >= 2) {
151 out << "Modulation Vector 2: " << setw(12) << setprecision(4) << lattice.getdh(1) << setw(12) << setprecision(4)
152 << lattice.getdk(1) << setw(12) << setprecision(4) << lattice.getdl(1) << " \n";
153
154 out << "Modulation Vector 2 error: " << setw(6) << setprecision(4) << lattice.getdherr(1) << setw(12)
155 << setprecision(4) << lattice.getdkerr(1) << setw(12) << setprecision(4) << lattice.getdlerr(1) << " \n";
156 }
157 if (ModDim == 3) {
158 out << "Modulation Vector 3: " << setw(12) << setprecision(4) << lattice.getdh(2) << setw(12) << setprecision(4)
159 << lattice.getdk(2) << setw(12) << setprecision(4) << lattice.getdl(2) << " \n";
160
161 out << "Modulation Vector 3 error: " << setw(6) << setprecision(4) << lattice.getdherr(2) << setw(12)
162 << setprecision(4) << lattice.getdkerr(2) << setw(12) << setprecision(4) << lattice.getdlerr(2) << " \n";
163 }
164 if (ModDim >= 1) {
165 out << "\n";
166 out << "Max Order: " << lattice.getMaxOrder() << " \n";
167 out << "Cross Terms: " << lattice.getCrossTerm() << " \n";
168 }
169
170 out << "\n";
171
172 if (ModDim == 0) {
173 out << "The above matrix is the Transpose of the UB Matrix. ";
174 out << "The UB matrix maps the column\n";
175 out << "vector (h,k,l ) to the column vector ";
176 out << "(q'x,q'y,q'z).\n";
177 out << "|Q'|=1/dspacing and its coordinates are a ";
178 out << "right-hand coordinate system where\n";
179 out << " x is the beam direction and z is vertically ";
180 out << "upward.(IPNS convention)\n";
181 } else {
182 out << "The above matrix is the Transpose of the UB Matrix and the "
183 "Transpose of ModUB. ";
184 out << "The UB matrix together with ModUB maps the column vector "
185 "(h,k,l,m,n,p) \n";
186 out << "to the column vector (q'x,q'y,q'z).\n";
187 out << "The columns of ModUB are the coordinates of modulation vectors "
188 "in Qlab. \n";
189 out << "|Q'|=1/dspacing and its coordinates are a ";
190 out << "right-hand coordinate system where";
191 out << " x is the beam direction and z is vertically ";
192 out << "upward.(IPNS convention)\n";
193 }
194
195 out.close();
196
197 } catch (exception &s) {
198 throw std::invalid_argument(s.what());
199 }
200}
201
202} // namespace Mantid::Crystal
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
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
@ Save
to specify a file to write to, the file may or may not exist
Definition: FileProperty.h:49
A property class for workspaces.
void init() override
Initialise the properties.
Definition: SaveIsawUB.cpp:29
void exec() override
Run the algorithm.
Definition: SaveIsawUB.cpp:77
double getErrorVolume(const Geometry::OrientedLattice &lattice)
Definition: SaveIsawUB.cpp:38
Class to implement UB matrix.
const Kernel::DblMatrix & getUB() const
Get the UB matrix.
const Kernel::DblMatrix & getModUB() const
double alpha() const
Get lattice parameter.
Definition: UnitCell.cpp:133
double erroralpha(const int angleunit=angDegrees) const
Get lattice parameter error.
Definition: UnitCell.cpp:225
double a(int nd) const
Get lattice parameter a1-a3 as function of index (0-2)
Definition: UnitCell.cpp:94
int getMaxOrder() const
Get max order.
Definition: UnitCell.cpp:596
double c() const
Get lattice parameter.
Definition: UnitCell.cpp:128
double getdh(int j) const
Get modulation vectors for satellites.
Definition: UnitCell.cpp:555
double getdk(int j) const
Get modulation vectors for satellites.
Definition: UnitCell.cpp:562
double volume() const
Volume of the direct unit-cell.
Definition: UnitCell.cpp:737
bool getCrossTerm() const
Get cross term boolean.
Definition: UnitCell.cpp:602
double getdl(int j) const
Get modulation vectors for satellites.
Definition: UnitCell.cpp:569
double errorgamma(const int angleunit=angDegrees) const
Get lattice parameter error.
Definition: UnitCell.cpp:252
double getdherr(int j) const
Get error of modulation vectors for satellites.
Definition: UnitCell.cpp:576
double errorbeta(const int angleunit=angDegrees) const
Get lattice parameter error.
Definition: UnitCell.cpp:239
double beta() const
Get lattice parameter.
Definition: UnitCell.cpp:138
double getdlerr(int j) const
Get error of modulation vectors for satellites.
Definition: UnitCell.cpp:590
double b() const
Get lattice parameter.
Definition: UnitCell.cpp:123
double getdkerr(int j) const
Get error of modulation vectors for satellites.
Definition: UnitCell.cpp:583
const Kernel::V3D getModVec(int j) const
Get modulation vectors for satellites.
Definition: UnitCell.cpp:535
double errorc() const
Get lattice parameter error.
Definition: UnitCell.cpp:218
double gamma() const
Get lattice parameter.
Definition: UnitCell.cpp:143
double errora() const
Get lattice parameter error.
Definition: UnitCell.cpp:208
double errorb() const
Get lattice parameter error.
Definition: UnitCell.cpp:213
Class for 3D vectors.
Definition: V3D.h:34
std::shared_ptr< Workspace > Workspace_sptr
shared pointer to Mantid::API::Workspace
Definition: Workspace_fwd.h:20
std::shared_ptr< ExperimentInfo > ExperimentInfo_sptr
Shared pointer to ExperimentInfo.
std::shared_ptr< MultipleExperimentInfos > MultipleExperimentInfos_sptr
STL namespace.
@ Input
An input workspace.
Definition: Property.h:53