Mantid
Loading...
Searching...
No Matches
OptimizeCrystalPlacement.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 +
7/*
8 *
9 * OptimizeCrystalPlacement.cpp
10 *
11 * Created on: Jan 26, 2013
12 * Author: ruth
13 */
15
16#include "MantidAPI/Run.h"
17#include "MantidAPI/Sample.h"
28
29#include <cstdarg>
30#include <utility>
31
32using namespace Mantid::API;
33using namespace Mantid::DataObjects;
34using namespace Mantid::Kernel;
37using namespace Mantid::Geometry;
38
39namespace Mantid::Crystal {
40
41DECLARE_ALGORITHM(OptimizeCrystalPlacement)
42
44public:
45 OrEnabledWhenProperties(std::string prop1Name, ePropertyCriterion prop1Crit, std::string prop1Value,
46 std::string prop2Name, ePropertyCriterion prop2Crit, std::string prop2Value)
47 : IPropertySettings(), propName1(std::move(prop1Name)), propName2(std::move(prop2Name)), Criteria1(prop1Crit),
48 Criteria2(prop2Crit), value1(std::move(prop1Value)), value2(std::move(prop2Value))
49
50 {
51 Prop1 = std::make_unique<Kernel::EnabledWhenProperty>(propName1, Criteria1, value1);
52 Prop2 = std::make_unique<Kernel::EnabledWhenProperty>(propName2, Criteria2, value2);
53 }
54 ~OrEnabledWhenProperties() override // responsible for deleting all supplied
55 // EnabledWhenProperites
56 = default;
57
58 IPropertySettings *clone() const override {
59 return new OrEnabledWhenProperties(propName1, Criteria1, value1, propName2, Criteria2, value2);
60 }
61
62 bool isEnabled(const IPropertyManager *algo) const override {
63 return Prop1->isEnabled(algo) && Prop2->isEnabled(algo);
64 }
65
66private:
67 std::string propName1, propName2;
69 std::string value1, value2;
70 std::unique_ptr<Kernel::EnabledWhenProperty> Prop1, Prop2;
71};
72
74 declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace>>("PeaksWorkspace", "", Direction::Input),
75 "Workspace of Peaks with UB loaded");
76 declareProperty(std::make_unique<ArrayProperty<int>>(std::string("KeepGoniometerFixedfor"), Direction::Input),
77 "List of run Numbers for which the goniometer settings will "
78 "NOT be changed");
79
80 declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace>>("ModifiedPeaksWorkspace", "", Direction::Output),
81 "Output Workspace of Peaks with optimized sample Orientations");
82
84 std::make_unique<WorkspaceProperty<ITableWorkspace>>("FitInfoTable", "FitInfoTable", Direction::Output),
85 "Workspace of Results");
86
87 declareProperty("AdjustSampleOffsets", false,
88 "If true sample offsets will be adjusted to give better "
89 "fits, otherwise they will be fixed as zero(def=true)");
90
91 declareProperty("OptimizeGoniometerTilt", false, "Set true if main error is due to a tilted Goniometer(def=false)");
92
93 declareProperty("Chi2overDoF", -1.0, "chi squared over dof", Direction::Output);
94 declareProperty("nPeaks", -1, "Number of Peaks Used", Direction::Output);
95 declareProperty("nParams", -1, "Number of Parameters fit", Direction::Output);
96 declareProperty("nIndexed", -1, "Number of new Peaks that WOULD be indexed at 'MaxIndexingError'", Direction::Output);
97
98 declareProperty("MaxAngularChange", 5.0, "Max offset in degrees from current settings(def=5)");
99
100 declareProperty("MaxIndexingError", 0.15,
101 "Use only peaks whose fractional "
102 "hkl values are below this "
103 "tolerance(def=0.15)");
104
105 declareProperty("MaxHKLPeaks2Use", -1.0,
106 "If less than 0 all peaks are used, "
107 "otherwise only peaks whose h,k, "
108 "and l values are below the level "
109 "are used(def=-1)");
110 declareProperty("MaxSamplePositionChangeMeters", .0005, "Maximum Change in Sample position in meters(def=.0005)");
111
112 setPropertyGroup("MaxAngularChange", "Tolerance settings");
113
114 setPropertyGroup("MaxSamplePositionChangeMeters", "Tolerance settings");
115 setPropertyGroup("MaxHKLPeaks2Use", "Tolerance settings");
116 setPropertyGroup("MaxIndexingError", "Tolerance settings");
117
118 setPropertySettings("MaxSamplePositionChangeMeters",
119 std::make_unique<EnabledWhenProperty>("AdjustSampleOffsets", Kernel::IS_EQUAL_TO, "1"));
120
121 setPropertySettings("KeepGoniometerFixedfor",
122 std::make_unique<OrEnabledWhenProperties>("AdjustSampleOffsets", Kernel::IS_EQUAL_TO, "0",
123 "OptimizeGoniometerTilt", Kernel::IS_EQUAL_TO, "0"));
124
125 declareProperty(std::make_unique<WorkspaceProperty<ITableWorkspace>>("OutputNormalisedCovarianceMatrixOptX",
126 "CovarianceInfo", Direction::Output),
127 "The name of the TableWorkspace in which to store the final "
128 "covariance matrix");
129}
130
140 PeaksWorkspace_sptr peaks = getProperty("PeaksWorkspace");
141 PeaksWorkspace_sptr outPeaks = getProperty("ModifiedPeaksWorkspace");
142
143 if (peaks != outPeaks) {
144 outPeaks = peaks->clone();
145 }
146
147 std::vector<int> NOoptimizeRuns = getProperty("KeepGoniometerFixedfor");
148
149 const DblMatrix X = peaks->sample().getOrientedLattice().getUB();
150 Matrix<double> UBinv(X);
151 UBinv.Invert();
152
153 //--------------------------------- Set up data for call to PeakHKLErrors
154 // fitting function ----------
155 // ---- Setting up workspace supplied to PeakHKLErrors
156 // ---------------
157 std::vector<int> RunNumList;
158 std::vector<V3D> ChiPhiOmega;
160
161 int nPeaksUsed = 0;
162 double HKLintOffsetMax = getProperty("MaxIndexingError");
163 double HKLMax = getProperty("MaxHKLPeaks2Use");
164 for (int i = 0; i < peaks->getNumberPeaks(); i++) {
165 IPeak &peak = peaks->getPeak(i);
166 int runNum = peak.getRunNumber();
167 auto it = RunNumList.begin();
168 for (; it != RunNumList.end() && *it != runNum; ++it) {
169 }
170
171 V3D hkl = UBinv * (peak.getQSampleFrame()) / (2.0 * M_PI);
172 bool use = IndexingUtils::ValidIndex(hkl, HKLintOffsetMax); // use this peak???
173
174 if (use && HKLMax > 0)
175 for (int k = 0; k < 3; k++) {
176 if (fabs(hkl[k]) > HKLMax)
177 use = false;
178 }
179
180 if (it == RunNumList.end() && use) // add to list of unique run numbers in workspace
181 {
182 RunNumList.emplace_back(runNum);
183
185 std::vector<double> phichiOmega = Gon.getEulerAngles("YZY");
186 ChiPhiOmega.emplace_back(phichiOmega[1], phichiOmega[2], phichiOmega[0]);
187 }
188
189 if (use) // add to lists for workspace
190 {
191 nPeaksUsed++;
192 xRef.emplace_back(static_cast<double>(i));
193 xRef.emplace_back(static_cast<double>(i));
194 xRef.emplace_back(static_cast<double>(i));
195 }
196 }
197
198 g_log.notice() << "Number initially indexed = " << nPeaksUsed << " at tolerance = " << HKLintOffsetMax << '\n';
199
200 if (nPeaksUsed < 1) {
201 g_log.error() << "Error in UB too large. 0 peaks indexed at " << HKLintOffsetMax << '\n';
202 throw std::invalid_argument("Error in UB too large. 0 peaks indexed ");
203 }
204
205 int N = 3 * nPeaksUsed; // peaks->getNumberPeaks();
206 auto mwkspc = createWorkspace<Workspace2D>(1, N, N);
207 mwkspc->setPoints(0, xRef);
208 mwkspc->setCounts(0, N, 0.0);
209 mwkspc->setCountStandardDeviations(0, N, 1.0);
210
211 std::string FuncArg = "name=PeakHKLErrors,PeakWorkspaceName=" + getPropertyValue("PeaksWorkspace") + "";
212
213 std::string OptRunNums;
214
215 //--------- Setting Function and Constraint argumens to PeakHKLErrors
216 //---------------
217 std::vector<std::string> ChRunNumList;
218 std::string predChar;
219 for (auto runNum : RunNumList) {
220 auto it1 = NOoptimizeRuns.begin();
221 for (; it1 != NOoptimizeRuns.end() && *it1 != runNum; ++it1) {
222 }
223
224 if (it1 == NOoptimizeRuns.end()) {
225 std::string runNumStr = std::to_string(runNum);
226 OptRunNums += predChar + runNumStr;
227 predChar = "/";
228 ChRunNumList.emplace_back(runNumStr);
229 }
230 }
231
232 bool omitRuns = (bool)getProperty("AdjustSampleOffsets") || (bool)getProperty("OptimizeGoniometerTilt");
233 if (omitRuns) {
234 NOoptimizeRuns = RunNumList;
235 OptRunNums = "";
236 std::string message = "No Goniometer Angles ";
237 if ((bool)getProperty("OptimizeGoniometerTilt"))
238 message += "relative to the tilted Goniometer ";
239 message += "will be 'changed'";
240 g_log.notice(message);
241 }
242 if (!OptRunNums.empty() && !omitRuns)
243 FuncArg += ",OptRuns=" + OptRunNums;
244
245 //------------- Add initial parameter values to FuncArg -----------
246
247 std::ostringstream oss(std::ostringstream::out);
248 oss.precision(3);
249 std::ostringstream oss1(std::ostringstream::out); // constraints
250 oss1.precision(3);
251
252 int nParams = 3;
253 double DegreeTol = getProperty("MaxAngularChange");
254 std::string startConstraint;
255
256 for (size_t i = 0; i < RunNumList.size(); i++) {
257 int runNum = RunNumList[i];
258
259 size_t k = 0;
260 for (; k < NOoptimizeRuns.size(); k++) {
261 if (NOoptimizeRuns[k] == runNum)
262 break;
263 }
264 if (k >= NOoptimizeRuns.size()) {
265 V3D chiphiomega = ChiPhiOmega[i];
266 oss << ",chi" << runNum << "=" << chiphiomega[0] << ",phi" << runNum << "=" << chiphiomega[1] << ",omega"
267 << runNum << "=" << chiphiomega[2];
268
269 oss1 << startConstraint << chiphiomega[0] - DegreeTol << "<chi" << runNum << "<" << chiphiomega[0] + DegreeTol;
270 oss1 << "," << chiphiomega[1] - DegreeTol << "<phi" << runNum << "<" << chiphiomega[1] + DegreeTol;
271
272 oss1 << "," << chiphiomega[2] - DegreeTol << "<omega" << runNum << "<" << chiphiomega[2] + DegreeTol;
273 startConstraint = ",";
274 nParams += 3;
275 }
276 }
277
278 // offset of previous sample position so should start at 0
279 V3D sampPos = V3D(0., 0., 0.);
280
281 oss << ",SampleXOffset=" << sampPos.X() << ",SampleYOffset=" << sampPos.Y() << ",SampleZOffset=" << sampPos.Z();
282 oss << ",GonRotx=0.0,GonRoty=0.0,GonRotz=0.0";
283
284 double maxSampshift = getProperty("MaxSamplePositionChangeMeters");
285 oss1 << startConstraint << sampPos.X() - maxSampshift << "<SampleXOffset<" << sampPos.X() + maxSampshift << ","
286 << sampPos.Y() - maxSampshift << "<SampleYOffset<" << sampPos.Y() + maxSampshift << ","
287 << sampPos.Z() - maxSampshift << "<SampleZOffset<" << sampPos.Z() + maxSampshift;
288
289 oss1 << "," << -DegreeTol << "<GonRotx<" << DegreeTol << "," << -DegreeTol << "<GonRoty<" << DegreeTol << ","
290 << -DegreeTol << "<GonRotz<" << DegreeTol;
291
292 FuncArg += oss.str();
293 std::string Constr = oss1.str();
294
295 g_log.debug() << "Function argument=" << FuncArg << '\n';
296 g_log.debug() << "Constraint argument=" << Constr << '\n';
297
298 //--------------------- set up Fit algorithm call-----------------
299
300 std::shared_ptr<Algorithm> fit_alg = createChildAlgorithm("Fit", .1, .93, true);
301
302 fit_alg->setProperty("Function", FuncArg);
303
304 fit_alg->setProperty("MaxIterations", 60);
305
306 fit_alg->setProperty("Constraints", Constr);
307
308 fit_alg->setProperty("InputWorkspace", mwkspc);
309
310 fit_alg->setProperty("CreateOutput", true);
311
312 std::string Ties;
313
314 if (!(bool)getProperty("AdjustSampleOffsets")) {
315 std::ostringstream oss3(std::ostringstream::out);
316 oss3.precision(3);
317
318 oss3 << "SampleXOffset=" << sampPos.X() << ",SampleYOffset=" << sampPos.Y() << ",SampleZOffset=" << sampPos.Z();
319 Ties = oss3.str();
320 }
321
322 if (!(bool)getProperty("OptimizeGoniometerTilt")) {
323 if (!Ties.empty())
324 Ties += ",";
325
326 Ties += "GonRotx=0.0,GonRoty=0.0,GonRotz=0.0";
327 }
328
329 if (!Ties.empty())
330 fit_alg->setProperty("Ties", Ties);
331
332 fit_alg->setProperty("Output", "out");
333
334 fit_alg->executeAsChildAlg();
335
336 //------------------------- Get/Report Results ------------------
337
338 double chisq = fit_alg->getProperty("OutputChi2overDoF");
339 g_log.notice() << "Fit finished. Status=" << (std::string)fit_alg->getProperty("OutputStatus") << '\n';
340
341 setProperty("Chi2overDoF", chisq);
342
343 setProperty("nPeaks", nPeaksUsed);
344 setProperty("nParams", nParams);
345
346 g_log.debug() << "Chi2overDof=" << chisq << " # Peaks used=" << nPeaksUsed << "# fitting parameters =" << nParams
347 << " dof=" << (nPeaksUsed - nParams) << '\n';
348
349 ITableWorkspace_sptr RRes = fit_alg->getProperty("OutputParameters");
350
351 double sigma = sqrt(chisq);
352
353 std::string OutputStatus = fit_alg->getProperty("OutputStatus");
354 g_log.notice() << "Output Status=" << OutputStatus << '\n';
355
356 //------------------ Fix up Covariance output --------------------
357 ITableWorkspace_sptr NormCov = fit_alg->getProperty("OutputNormalisedCovarianceMatrix");
358 setProperty("OutputNormalisedCovarianceMatrixOptX",
359 NormCov); // What if 2 instances are run
360
361 if (chisq < 0 || chisq != chisq)
362 sigma = -1;
363
364 //------------- Fix up Result workspace values ----------------------------
365 std::map<std::string, double> Results;
366 for (int prm = 0; prm < static_cast<int>(RRes->rowCount()); ++prm) {
367 std::string namee = RRes->getRef<std::string>("Name", prm);
368
369 std::string start = namee.substr(0, 3);
370 if (start != "chi" && start != "phi" && start != "ome" && start != "Sam" && start != "Gon")
371 continue;
372
373 double value = RRes->getRef<double>("Value", prm);
374 Results[namee] = value;
375
376 // Set sigma==1 in optimization. A better estimate is sqrt(Chi2overDoF)
377 double v = sigma * RRes->getRef<double>("Error", prm);
378 RRes->getRef<double>("Error", prm) = v;
379 }
380
381 //-----------Fix up Resultant workspace return info -------------------
382
383 setProperty("FitInfoTable", RRes);
384
385 //----------- update instrument -------------------------
386 V3D newSampPos(Results["SampleXOffset"], Results["SampleYOffset"], Results["SampleZOffset"]);
387
388 auto &componentInfo = outPeaks->mutableComponentInfo();
389 CalibrationHelpers::adjustUpSampleAndSourcePositions(componentInfo.l1(), newSampPos, componentInfo);
390
391 Matrix<double> GonTilt = PeakHKLErrors::RotationMatrixAboutRegAxis(Results["GonRotx"], 'x') *
392 PeakHKLErrors::RotationMatrixAboutRegAxis(Results["GonRoty"], 'y') *
393 PeakHKLErrors::RotationMatrixAboutRegAxis(Results["GonRotz"], 'z');
394
395 int prevRunNum = -1;
396 std::map<int, Matrix<double>> MapRunNum2GonMat;
397 std::string OptRun2 = "/" + OptRunNums + "/";
398
399 int nIndexed = 0;
400 UBinv = outPeaks->sample().getOrientedLattice().getUB();
401 UBinv.Invert();
402 UBinv /= (2 * M_PI);
403 for (int i = 0; i < outPeaks->getNumberPeaks(); ++i) {
404 auto &peak = outPeaks->getPeak(i);
405 peak.setSamplePos(peak.getSamplePos() + newSampPos);
406 int RunNum = peak.getRunNumber();
407 std::string RunNumStr = std::to_string(RunNum);
408 Matrix<double> GonMatrix;
409 if (RunNum == prevRunNum || MapRunNum2GonMat.find(RunNum) != MapRunNum2GonMat.end())
410 GonMatrix = MapRunNum2GonMat[RunNum];
411 else if (OptRun2.find("/" + RunNumStr + "/") < OptRun2.size() - 2) {
412
413 double chi = Results["chi" + RunNumStr];
414 double phi = Results["phi" + RunNumStr];
415 double omega = Results["omega" + RunNumStr];
416
418 uniGonio.makeUniversalGoniometer();
419
420 uniGonio.setRotationAngle("phi", phi);
421 uniGonio.setRotationAngle("chi", chi);
422 uniGonio.setRotationAngle("omega", omega);
423
424 GonMatrix = GonTilt * uniGonio.getR();
425 MapRunNum2GonMat[RunNum] = GonMatrix;
426 } else {
427 GonMatrix = GonTilt * peak.getGoniometerMatrix();
428 MapRunNum2GonMat[RunNum] = GonMatrix;
429 }
430
431 peak.setGoniometerMatrix(GonMatrix);
432 V3D hkl = UBinv * peak.getQSampleFrame();
433 if (Geometry::IndexingUtils::ValidIndex(hkl, HKLintOffsetMax))
434 nIndexed++;
435
436 prevRunNum = RunNum;
437 }
438
439 if (MapRunNum2GonMat.size() == 1) // Only one RunNumber in this PeaksWorkspace
440 {
441 Matrix<double> GonMatrix = MapRunNum2GonMat[outPeaks->getPeak(0).getRunNumber()];
442 Geometry::Goniometer Gon(GonMatrix);
443 outPeaks->mutableRun().setGoniometer(Gon, false);
444 }
445
446 setProperty("ModifiedPeaksWorkspace", outPeaks);
447 setProperty("nIndexed", nIndexed);
448 g_log.notice() << "Number indexed after optimization= " << nIndexed << " at tolerance = " << HKLintOffsetMax << '\n';
449
450} // exec
451
452} // namespace Mantid::Crystal
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double value
The value of the point.
Definition: FitMW.cpp:51
double sigma
Definition: GetAllEi.cpp:156
#define fabs(x)
Definition: Matrix.cpp:22
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
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
Definition: Algorithm.cpp:2026
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
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
A property class for workspaces.
void init() override
Virtual method - must be overridden by concrete algorithm.
bool isEnabled(const IPropertyManager *algo) const override
Is the property to be shown as "enabled" in the GUI.
std::unique_ptr< Kernel::EnabledWhenProperty > Prop1
OrEnabledWhenProperties(std::string prop1Name, ePropertyCriterion prop1Crit, std::string prop1Value, std::string prop2Name, ePropertyCriterion prop2Crit, std::string prop2Value)
static Kernel::Matrix< double > RotationMatrixAboutRegAxis(double theta, char axis)
Returns the matrix corresponding to a rotation of theta(degrees) around axis.
Class to represent a particular goniometer setting, which is described by the rotation matrix.
Definition: Goniometer.h:55
void setRotationAngle(const std::string &name, double value)
Set rotation angle for an axis using motor name.
Definition: Goniometer.cpp:161
std::vector< double > getEulerAngles(const std::string &convention="YZX")
Return Euler angles acording to a convention.
Definition: Goniometer.cpp:282
void makeUniversalGoniometer()
Make a default universal goniometer with phi,chi,omega angles according to SNS convention.
Definition: Goniometer.cpp:271
const Kernel::DblMatrix & getR() const
Return global rotation matrix.
Definition: Goniometer.cpp:80
Structure describing a single-crystal peak.
Definition: IPeak.h:26
virtual Mantid::Kernel::Matrix< double > getGoniometerMatrix() const =0
virtual int getRunNumber() const =0
virtual Mantid::Kernel::V3D getQSampleFrame() const =0
This class contains static utility methods for indexing peaks and finding the UB matrix.
Definition: IndexingUtils.h:32
static bool ValidIndex(const Kernel::V3D &hkl, double tolerance)
Check is hkl is within tolerance of integer (h,k,l) non-zero values.
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
Interface to PropertyManager.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void setPropertySettings(const std::string &name, std::unique_ptr< IPropertySettings > settings)
void setPropertyGroup(const std::string &name, const std::string &group)
Set the group for a given property.
Interface for modifiers to Property's that specify if they should be enabled or visible in a GUI.
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
T Invert()
LU inversion routine.
Definition: Matrix.cpp:924
Class for 3D vectors.
Definition: V3D.h:34
constexpr double X() const noexcept
Get x.
Definition: V3D.h:232
constexpr double Y() const noexcept
Get y.
Definition: V3D.h:233
constexpr double Z() const noexcept
Get z.
Definition: V3D.h:234
std::shared_ptr< ITableWorkspace > ITableWorkspace_sptr
shared pointer to Mantid::API::ITableWorkspace
MANTID_CRYSTAL_DLL void adjustUpSampleAndSourcePositions(double const L0, const Kernel::V3D &newSampPos, Geometry::ComponentInfo &componentInfo)
Updates the ComponentInfo for the workspace containing newInstrument to reflect the position of the s...
std::shared_ptr< PeaksWorkspace > PeaksWorkspace_sptr
Typedef for a shared pointer to a peaks workspace.
std::shared_ptr< const Instrument > Instrument_const_sptr
Shared pointer to an const instrument object.
ePropertyCriterion
Enum for use in EnabledWhenProperty.
std::vector< double > MantidVec
typedef for the data storage used in Mantid matrix workspaces
Definition: cow_ptr.h:172
STL namespace.
std::string to_string(const wide_integer< Bits, Signed > &n)
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54