Skip to content
Snippets Groups Projects
Commit a5008c56 authored by Savici, Andrei T.'s avatar Savici, Andrei T. Committed by GitHub
Browse files

Merge pull request #18107 from mantidproject/18039_findublatticeparams_new_options 

Update FindUBUsingLatticeParameters for small numbers of peaks
parents 54730946 925940bd
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 224 additions and 109 deletions
Framework/Crystal/inc/MantidCrystal/FindUBUsingLatticeParameters.h
+ 1
− 1
View file @ a5008c56
#ifndef MANTID_CRYSTAL_FIND_UB_USING_LATTICE_PARAMETERS_H_ #ifndef MANTID_CRYSTAL_FIND_UB_USING_LATTICE_PARAMETERS_H_
#define MANTID_CRYSTAL_FIND_UB_USING_LATTICE_PARAMETERS_H_ #define MANTID_CRYSTAL_FIND_UB_USING_LATTICE_PARAMETERS_H_
#include "MantidKernel/System.h"
#include "MantidAPI/Algorithm.h" #include "MantidAPI/Algorithm.h"
#include "MantidKernel/System.h"
namespace Mantid { namespace Mantid {
namespace Crystal { namespace Crystal {
......
Framework/Crystal/src/FindUBUsingLatticeParameters.cpp
+ 17
− 31
View file @ a5008c56
...@@ -46,6 +46,8 @@ void FindUBUsingLatticeParameters::init() { ...@@ -46,6 +46,8 @@ void FindUBUsingLatticeParameters::init() {
"Number of Peaks to Use on First Pass(15)"); "Number of Peaks to Use on First Pass(15)");
this->declareProperty("Tolerance", 0.15, mustBePositive, this->declareProperty("Tolerance", 0.15, mustBePositive,
"Indexing Tolerance (0.15)"); "Indexing Tolerance (0.15)");
this->declareProperty("FixParameters", false,
"Do not optimise the UB after finding the orientation");
} }
/** Execute the algorithm. /** Execute the algorithm.
...@@ -59,6 +61,7 @@ void FindUBUsingLatticeParameters::exec() { ...@@ -59,6 +61,7 @@ void FindUBUsingLatticeParameters::exec() {
double gamma = this->getProperty("gamma"); double gamma = this->getProperty("gamma");
int num_initial = this->getProperty("NumInitial"); int num_initial = this->getProperty("NumInitial");
double tolerance = this->getProperty("Tolerance"); double tolerance = this->getProperty("Tolerance");
auto fixAll = this->getProperty("FixParameters");
int base_index = -1; // these "could" be properties if need be int base_index = -1; // these "could" be properties if need be
double degrees_per_step = 1.5; double degrees_per_step = 1.5;
...@@ -77,12 +80,13 @@ void FindUBUsingLatticeParameters::exec() { ...@@ -77,12 +80,13 @@ void FindUBUsingLatticeParameters::exec() {
q_vectors.push_back(peaks[i].getQSampleFrame()); q_vectors.push_back(peaks[i].getQSampleFrame());
Matrix<double> UB(3, 3, false); Matrix<double> UB(3, 3, false);
double error = IndexingUtils::Find_UB(UB, q_vectors, a, b, c, alpha, beta, OrientedLattice lattice(a, b, c, alpha, beta, gamma);
gamma, tolerance, base_index, double error =
num_initial, degrees_per_step); IndexingUtils::Find_UB(UB, q_vectors, lattice, tolerance, base_index,
num_initial, degrees_per_step, fixAll);
std::cout << "Error = " << error << '\n'; g_log.notice() << "Error = " << error << '\n';
std::cout << "UB = " << UB << '\n'; g_log.notice() << "UB = " << UB << '\n';
if (!IndexingUtils::CheckUB(UB)) // UB not found correctly if (!IndexingUtils::CheckUB(UB)) // UB not found correctly
{ {
...@@ -91,18 +95,6 @@ void FindUBUsingLatticeParameters::exec() { ...@@ -91,18 +95,6 @@ void FindUBUsingLatticeParameters::exec() {
g_log.notice(std::string("UB NOT SAVED.")); g_log.notice(std::string("UB NOT SAVED."));
} else // tell user how many would be indexed } else // tell user how many would be indexed
{ // and save the UB in the sample { // and save the UB in the sample
std::vector<double> sigabc(7);
std::vector<V3D> miller_ind;
std::vector<V3D> indexed_qs;
double fit_error;
miller_ind.reserve(q_vectors.size());
indexed_qs.reserve(q_vectors.size());
IndexingUtils::GetIndexedPeaks(UB, q_vectors, tolerance, miller_ind,
indexed_qs, fit_error);
IndexingUtils::Optimize_UB(UB, miller_ind, indexed_qs, sigabc);
char logInfo[200]; char logInfo[200];
int num_indexed = IndexingUtils::NumberIndexed(UB, q_vectors, tolerance); int num_indexed = IndexingUtils::NumberIndexed(UB, q_vectors, tolerance);
sprintf(logInfo, sprintf(logInfo,
...@@ -112,27 +104,21 @@ void FindUBUsingLatticeParameters::exec() { ...@@ -112,27 +104,21 @@ void FindUBUsingLatticeParameters::exec() {
num_indexed, n_peaks, tolerance); num_indexed, n_peaks, tolerance);
g_log.notice(std::string(logInfo)); g_log.notice(std::string(logInfo));
OrientedLattice o_lattice; double calc_a = lattice.a();
o_lattice.setUB(UB); double calc_b = lattice.b();
o_lattice.setError(sigabc[0], sigabc[1], sigabc[2], sigabc[3], sigabc[4], double calc_c = lattice.c();
sigabc[5]); double calc_alpha = lattice.alpha();
double calc_beta = lattice.beta();
o_lattice.setUB(UB); double calc_gamma = lattice.gamma();
double calc_a = o_lattice.a();
double calc_b = o_lattice.b();
double calc_c = o_lattice.c();
double calc_alpha = o_lattice.alpha();
double calc_beta = o_lattice.beta();
double calc_gamma = o_lattice.gamma();
// Show the modified lattice parameters // Show the modified lattice parameters
g_log.notice() << o_lattice << "\n"; g_log.notice() << lattice << "\n";
sprintf(logInfo, std::string("Lattice Parameters (Refined - Input): %11.6f " sprintf(logInfo, std::string("Lattice Parameters (Refined - Input): %11.6f "
"%11.6f %11.6f %11.6f %11.6f %11.6f").c_str(), "%11.6f %11.6f %11.6f %11.6f %11.6f").c_str(),
calc_a - a, calc_b - b, calc_c - c, calc_alpha - alpha, calc_a - a, calc_b - b, calc_c - c, calc_alpha - alpha,
calc_beta - beta, calc_gamma - gamma); calc_beta - beta, calc_gamma - gamma);
g_log.notice(std::string(logInfo)); g_log.notice(std::string(logInfo));
ws->mutableSample().setOrientedLattice(&o_lattice); ws->mutableSample().setOrientedLattice(&lattice);
} }
} }
......
Framework/Crystal/test/FindUBUsingLatticeParametersTest.h
+ 152
− 22
View file @ a5008c56
#ifndef MANTID_CRYSTAL_FIND_UB_USING_LATTICE_PARAMETERS_TEST_H_ #ifndef MANTID_CRYSTAL_FIND_UB_USING_LATTICE_PARAMETERS_TEST_H_
#define MANTID_CRYSTAL_FIND_UB_USING_LATTICE_PARAMETERS_TEST_H_ #define MANTID_CRYSTAL_FIND_UB_USING_LATTICE_PARAMETERS_TEST_H_
#include <cxxtest/TestSuite.h>
#include "MantidKernel/Timer.h"
#include "MantidKernel/System.h"
#include "MantidAPI/Sample.h" #include "MantidAPI/Sample.h"
#include "MantidDataHandling/DeleteTableRows.h"
#include "MantidKernel/System.h"
#include "MantidKernel/Timer.h"
#include <cxxtest/TestSuite.h>
#include "MantidCrystal/FindUBUsingLatticeParameters.h" #include "MantidCrystal/FindUBUsingLatticeParameters.h"
#include "MantidCrystal/LoadIsawPeaks.h" #include "MantidCrystal/LoadIsawPeaks.h"
...@@ -27,28 +28,13 @@ public: ...@@ -27,28 +28,13 @@ public:
} }
void test_exec() { void test_exec() {
// Name of the output workspace. auto ws = loadPeaksWorkspace();
std::string WSName("peaks");
LoadIsawPeaks loader;
TS_ASSERT_THROWS_NOTHING(loader.initialize());
TS_ASSERT(loader.isInitialized());
loader.setPropertyValue("Filename", "TOPAZ_3007.peaks");
loader.setPropertyValue("OutputWorkspace", WSName);
TS_ASSERT(loader.execute());
TS_ASSERT(loader.isExecuted());
PeaksWorkspace_sptr ws;
TS_ASSERT_THROWS_NOTHING(
ws = boost::dynamic_pointer_cast<PeaksWorkspace>(
AnalysisDataService::Instance().retrieve(WSName)));
TS_ASSERT(ws);
if (!ws)
return;
FindUBUsingLatticeParameters alg; FindUBUsingLatticeParameters alg;
TS_ASSERT_THROWS_NOTHING(alg.initialize()) TS_ASSERT_THROWS_NOTHING(alg.initialize())
TS_ASSERT(alg.isInitialized()) TS_ASSERT(alg.isInitialized())
TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("PeaksWorkspace", WSName)); TS_ASSERT_THROWS_NOTHING(
alg.setPropertyValue("PeaksWorkspace", ws->name()));
TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("a", "14.131")); TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("a", "14.131"));
TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("b", "19.247")); TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("b", "19.247"));
TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("c", "8.606")); TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("c", "8.606"));
...@@ -75,8 +61,152 @@ public: ...@@ -75,8 +61,152 @@ public:
TS_ASSERT_DELTA(correct_UB[i], UB_calculated[i], 5e-4); TS_ASSERT_DELTA(correct_UB[i], UB_calculated[i], 5e-4);
} }
TS_ASSERT_DELTA(latt.a(), 14.131, 5e-4);
TS_ASSERT_DELTA(latt.b(), 19.247, 5e-4);
TS_ASSERT_DELTA(latt.c(), 8.606, 5e-4);
TS_ASSERT_DELTA(latt.alpha(), 90.0, 5e-1);
TS_ASSERT_DELTA(latt.beta(), 105.071, 5e-1);
TS_ASSERT_DELTA(latt.gamma(), 90.0, 5e-1);
// Check errors
TS_ASSERT_DELTA(latt.errora(), 0.0134, 5e-4);
TS_ASSERT_DELTA(latt.errorb(), 0.0243, 5e-4);
TS_ASSERT_DELTA(latt.errorc(), 0.0101, 5e-4);
TS_ASSERT_DELTA(latt.erroralpha(), 0.0994, 5e-4);
TS_ASSERT_DELTA(latt.errorbeta(), 0.0773, 5e-4);
TS_ASSERT_DELTA(latt.errorgamma(), 0.0906, 5e-4);
// Remove workspace from the data service.
AnalysisDataService::Instance().remove(ws->name());
}
void test_fixAll() {
auto ws = loadPeaksWorkspace();
FindUBUsingLatticeParameters alg;
TS_ASSERT_THROWS_NOTHING(alg.initialize())
TS_ASSERT(alg.isInitialized())
TS_ASSERT_THROWS_NOTHING(
alg.setPropertyValue("PeaksWorkspace", ws->name()));
TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("a", "14.131"));
TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("b", "19.247"));
TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("c", "8.606"));
TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("alpha", "90.0"));
TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("beta", "105.071"));
TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("gamma", "90.0"));
TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("NumInitial", "15"));
TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("Tolerance", "0.12"));
TS_ASSERT_THROWS_NOTHING(alg.setProperty("FixParameters", true));
TS_ASSERT_THROWS_NOTHING(alg.execute());
TS_ASSERT(alg.isExecuted());
// Check that we set an oriented lattice
TS_ASSERT(ws->mutableSample().hasOrientedLattice());
// Check that the UB matrix is the same as in TOPAZ_3007.mat
OrientedLattice latt = ws->mutableSample().getOrientedLattice();
double correct_UB[] = {0.04542050, 0.040619900, 0.0127661,
-0.00198382, -0.00264404, -0.1165450,
-0.05749760, 0.03223800, -0.0257623};
std::vector<double> UB_calculated = latt.getUB().getVector();
for (size_t i = 0; i < 9; i++) {
TS_ASSERT_DELTA(correct_UB[i], UB_calculated[i], 5e-4);
}
TS_ASSERT_DELTA(latt.a(), 14.131, 5e-10);
TS_ASSERT_DELTA(latt.b(), 19.247, 5e-10);
TS_ASSERT_DELTA(latt.c(), 8.606, 5e-10);
TS_ASSERT_DELTA(latt.alpha(), 90.0, 5e-10);
TS_ASSERT_DELTA(latt.beta(), 105.071, 5e-10);
TS_ASSERT_DELTA(latt.gamma(), 90.0, 5e-10);
// Remove workspace from the data service. // Remove workspace from the data service.
AnalysisDataService::Instance().remove(WSName); AnalysisDataService::Instance().remove(ws->name());
}
void test_smallNumberOfPeaks() {
// Use a tiny set of 3 peaks - the minimum required
/// to successfully find a UB matrix this checks the case that we still
/// get a UB (although perhaps not a very good one).
auto ws = loadPeaksWorkspace();
std::vector<size_t> rows;
for (size_t i = 3; i < ws->rowCount(); ++i) {
rows.push_back(i);
}
DataHandling::DeleteTableRows removeRowAlg;
removeRowAlg.initialize();
removeRowAlg.setPropertyValue("TableWorkspace", ws->name());
removeRowAlg.setProperty("Rows", rows);
removeRowAlg.execute();
FindUBUsingLatticeParameters alg;
TS_ASSERT_THROWS_NOTHING(alg.initialize())
TS_ASSERT(alg.isInitialized())
TS_ASSERT_THROWS_NOTHING(
alg.setPropertyValue("PeaksWorkspace", ws->name()));
TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("a", "14.131"));
TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("b", "19.247"));
TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("c", "8.606"));
TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("alpha", "90.0"));
TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("beta", "105.071"));
TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("gamma", "90.0"));
TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("NumInitial", "15"));
TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("Tolerance", "0.12"));
// TS_ASSERT_THROWS_NOTHING(alg.setProperty("FixAll", true));
TS_ASSERT_THROWS_NOTHING(alg.execute());
TS_ASSERT(alg.isExecuted());
// Check that we set an oriented lattice
TS_ASSERT(ws->mutableSample().hasOrientedLattice());
// Check that the UB matrix is the same as in TOPAZ_3007.mat
OrientedLattice latt = ws->mutableSample().getOrientedLattice();
double correct_UB[] = {0.0450, 0.0407, 0.0127, -0.0008, -0.0044,
-0.1158, -0.0584, 0.0307, -0.0242};
std::vector<double> UB_calculated = latt.getUB().getVector();
for (size_t i = 0; i < 9; i++) {
TS_ASSERT_DELTA(correct_UB[i], UB_calculated[i], 5e-4);
}
TS_ASSERT_DELTA(latt.a(), 13.9520, 5e-4);
TS_ASSERT_DELTA(latt.b(), 19.5145, 5e-4);
TS_ASSERT_DELTA(latt.c(), 8.6566, 5e-4);
TS_ASSERT_DELTA(latt.alpha(), 92.6267, 5e-4);
TS_ASSERT_DELTA(latt.beta(), 103.7440, 5e-4);
TS_ASSERT_DELTA(latt.gamma(), 90.0272, 5e-4);
// Remove workspace from the data service.
AnalysisDataService::Instance().remove(ws->name());
}
private:
/*
* Load a peaks workspace to use as input data
*/
PeaksWorkspace_sptr loadPeaksWorkspace() const {
std::string WSName("peaks");
LoadIsawPeaks loader;
TS_ASSERT_THROWS_NOTHING(loader.initialize());
TS_ASSERT(loader.isInitialized());
loader.setPropertyValue("Filename", "TOPAZ_3007.peaks");
loader.setPropertyValue("OutputWorkspace", WSName);
TS_ASSERT(loader.execute());
TS_ASSERT(loader.isExecuted());
PeaksWorkspace_sptr ws;
TS_ASSERT_THROWS_NOTHING(
ws = boost::dynamic_pointer_cast<PeaksWorkspace>(
AnalysisDataService::Instance().retrieve(WSName)));
TS_ASSERT(ws);
return ws;
} }
}; };
......
Framework/Geometry/inc/MantidGeometry/Crystal/IndexingUtils.h
+ 9
− 9
View file @ a5008c56
...@@ -6,10 +6,11 @@ ...@@ -6,10 +6,11 @@
//---------------------------------------------------------------------- //----------------------------------------------------------------------
// Includes // Includes
//---------------------------------------------------------------------- //----------------------------------------------------------------------
#include "MantidGeometry/Crystal/OrientedLattice.h"
#include "MantidGeometry/DllConfig.h" #include "MantidGeometry/DllConfig.h"
#include "MantidKernel/V3D.h"
#include "MantidKernel/Matrix.h"
#include "MantidKernel/Logger.h" #include "MantidKernel/Logger.h"
#include "MantidKernel/Matrix.h"
#include "MantidKernel/V3D.h"
namespace Mantid { namespace Mantid {
namespace Geometry { namespace Geometry {
...@@ -51,10 +52,10 @@ public: ...@@ -51,10 +52,10 @@ public:
/// Find the UB matrix that most nearly indexes the specified qxyz values /// Find the UB matrix that most nearly indexes the specified qxyz values
/// given the lattice parameters /// given the lattice parameters
static double Find_UB(Kernel::DblMatrix &UB, static double Find_UB(Kernel::DblMatrix &UB,
const std::vector<Kernel::V3D> &q_vectors, double a, const std::vector<Kernel::V3D> &q_vectors,
double b, double c, double alpha, double beta, OrientedLattice &lattice, double required_tolerance,
double gamma, double required_tolerance, int base_index, int base_index, size_t num_initial,
size_t num_initial, double degrees_per_step); double degrees_per_step, bool fixAll = false);
/// Find the UB matrix that most nearly indexes the specified qxyz values /// Find the UB matrix that most nearly indexes the specified qxyz values
/// given the range of possible real space unit cell edge lengths. /// given the range of possible real space unit cell edge lengths.
...@@ -88,9 +89,8 @@ public: ...@@ -88,9 +89,8 @@ public:
/// Scan rotations to find UB that indexes peaks given lattice parameters /// Scan rotations to find UB that indexes peaks given lattice parameters
static double ScanFor_UB(Kernel::DblMatrix &UB, static double ScanFor_UB(Kernel::DblMatrix &UB,
const std::vector<Kernel::V3D> &q_vectors, double a, const std::vector<Kernel::V3D> &q_vectors,
double b, double c, double alpha, double beta, const UnitCell &lattice, double degrees_per_step,
double gamma, double degrees_per_step,
double required_tolerance); double required_tolerance);
/// Get list of possible directions and lengths for real space unit cell /// Get list of possible directions and lengths for real space unit cell
......
Framework/Geometry/src/Crystal/IndexingUtils.cpp
+ 36
− 28
View file @ a5008c56
#include "MantidGeometry/Crystal/IndexingUtils.h" #include "MantidGeometry/Crystal/IndexingUtils.h"
#include "MantidGeometry/Crystal/NiggliCell.h" #include "MantidGeometry/Crystal/NiggliCell.h"
#include "MantidGeometry/Crystal/OrientedLattice.h"
#include "MantidKernel/Quat.h" #include "MantidKernel/Quat.h"
#include <algorithm> #include <algorithm>
#include <cmath> #include <cmath>
...@@ -53,12 +52,9 @@ const constexpr double RAD_TO_DEG = 180. / M_PI; ...@@ -53,12 +52,9 @@ const constexpr double RAD_TO_DEG = 180. / M_PI;
@param q_vectors std::vector of V3D objects that contains the @param q_vectors std::vector of V3D objects that contains the
list of q_vectors that are to be indexed list of q_vectors that are to be indexed
NOTE: There must be at least 2 q_vectors. NOTE: There must be at least 2 q_vectors.
@param a First unit cell edge length in Angstroms. @param lattice The orientated lattice with the lattice
@param b Second unit cell edge length in Angstroms. parameters a,b,c and alpha, beta, gamma. The found
@param c Third unit cell edge length in Angstroms. UB and errors will be set on this lattice.
@param alpha First unit cell angle in degrees.
@param beta second unit cell angle in degrees.
@param gamma third unit cell angle in degrees.
@param required_tolerance The maximum allowed deviation of Miller indices @param required_tolerance The maximum allowed deviation of Miller indices
from integer values for a peak to be indexed. from integer values for a peak to be indexed.
@param base_index The sequence number of the peak that should @param base_index The sequence number of the peak that should
...@@ -77,6 +73,8 @@ const constexpr double RAD_TO_DEG = 180. / M_PI; ...@@ -77,6 +73,8 @@ const constexpr double RAD_TO_DEG = 180. / M_PI;
used to scan for an initial orientation matrix. used to scan for an initial orientation matrix.
@param degrees_per_step The number of degrees between different @param degrees_per_step The number of degrees between different
orientations used during the initial scan. orientations used during the initial scan.
@param fixAll Fix the lattice parameters and do not optimise
the UB matrix.
@return This will return the sum of the squares of the residual errors. @return This will return the sum of the squares of the residual errors.
...@@ -87,10 +85,10 @@ const constexpr double RAD_TO_DEG = 180. / M_PI; ...@@ -87,10 +85,10 @@ const constexpr double RAD_TO_DEG = 180. / M_PI;
is <= 0. is <= 0.
*/ */
double IndexingUtils::Find_UB(DblMatrix &UB, const std::vector<V3D> &q_vectors, double IndexingUtils::Find_UB(DblMatrix &UB, const std::vector<V3D> &q_vectors,
double a, double b, double c, double alpha, OrientedLattice &lattice,
double beta, double gamma,
double required_tolerance, int base_index, double required_tolerance, int base_index,
size_t num_initial, double degrees_per_step) { size_t num_initial, double degrees_per_step,
bool fixAll) {
if (UB.numRows() != 3 || UB.numCols() != 3) { if (UB.numRows() != 3 || UB.numCols() != 3) {
throw std::invalid_argument("Find_UB(): UB matrix NULL or not 3X3"); throw std::invalid_argument("Find_UB(): UB matrix NULL or not 3X3");
} }
...@@ -156,19 +154,19 @@ double IndexingUtils::Find_UB(DblMatrix &UB, const std::vector<V3D> &q_vectors, ...@@ -156,19 +154,19 @@ double IndexingUtils::Find_UB(DblMatrix &UB, const std::vector<V3D> &q_vectors,
for (size_t i = 0; i < num_initial; i++) for (size_t i = 0; i < num_initial; i++)
some_qs.push_back(sorted_qs[i]); some_qs.push_back(sorted_qs[i]);
ScanFor_UB(UB, some_qs, a, b, c, alpha, beta, gamma, degrees_per_step, ScanFor_UB(UB, some_qs, lattice, degrees_per_step, required_tolerance);
required_tolerance);
double fit_error = 0; double fit_error = 0;
std::vector<V3D> miller_ind; std::vector<V3D> miller_ind;
std::vector<V3D> indexed_qs; std::vector<V3D> indexed_qs;
miller_ind.reserve(q_vectors.size()); miller_ind.reserve(q_vectors.size());
indexed_qs.reserve(q_vectors.size()); indexed_qs.reserve(q_vectors.size());
// now gradually bring in the remaining // now gradually bring in the remaining
// peaks and re-optimize the UB to index // peaks and re-optimize the UB to index
// them as well // them as well
size_t count = 0; size_t count = 0;
while (num_initial < sorted_qs.size()) { while (!fixAll && num_initial < sorted_qs.size()) {
count++; count++;
num_initial = std::lround(1.5 * static_cast<double>(num_initial + 3)); num_initial = std::lround(1.5 * static_cast<double>(num_initial + 3));
// add 3, in case we started with // add 3, in case we started with
...@@ -190,13 +188,15 @@ double IndexingUtils::Find_UB(DblMatrix &UB, const std::vector<V3D> &q_vectors, ...@@ -190,13 +188,15 @@ double IndexingUtils::Find_UB(DblMatrix &UB, const std::vector<V3D> &q_vectors,
} }
} }
if (q_vectors.size() >= 3) // try one last refinement using all peaks std::vector<double> sigabc(7);
if (!fixAll &&
q_vectors.size() >= 3) // try one last refinement using all peaks
{ {
try { try {
GetIndexedPeaks(UB, q_vectors, required_tolerance, miller_ind, indexed_qs, GetIndexedPeaks(UB, q_vectors, required_tolerance, miller_ind, indexed_qs,
fit_error); fit_error);
Matrix<double> temp_UB(3, 3, false); Matrix<double> temp_UB = UB;
fit_error = Optimize_UB(temp_UB, miller_ind, indexed_qs); fit_error = Optimize_UB(temp_UB, miller_ind, indexed_qs, sigabc);
UB = temp_UB; UB = temp_UB;
} catch (...) { } catch (...) {
// failed to improve UB using these peaks, so just return the current UB // failed to improve UB using these peaks, so just return the current UB
...@@ -207,6 +207,10 @@ double IndexingUtils::Find_UB(DblMatrix &UB, const std::vector<V3D> &q_vectors, ...@@ -207,6 +207,10 @@ double IndexingUtils::Find_UB(DblMatrix &UB, const std::vector<V3D> &q_vectors,
// HKL space. // HKL space.
GetIndexedPeaks(UB, q_vectors, required_tolerance, miller_ind, indexed_qs, GetIndexedPeaks(UB, q_vectors, required_tolerance, miller_ind, indexed_qs,
fit_error); fit_error);
// set the error on the lattice parameters
lattice.setUB(UB);
lattice.setError(sigabc[0], sigabc[1], sigabc[2], sigabc[3], sigabc[4],
sigabc[5]);
return fit_error; return fit_error;
} }
...@@ -583,7 +587,9 @@ double IndexingUtils::Optimize_UB(DblMatrix &UB, ...@@ -583,7 +587,9 @@ double IndexingUtils::Optimize_UB(DblMatrix &UB,
const std::vector<V3D> &hkl_vectors, const std::vector<V3D> &hkl_vectors,
const std::vector<V3D> &q_vectors, const std::vector<V3D> &q_vectors,
std::vector<double> &sigabc) { std::vector<double> &sigabc) {
double result = Optimize_UB(UB, hkl_vectors, q_vectors); double result = 0;
result = Optimize_UB(UB, hkl_vectors, q_vectors);
if (sigabc.size() < 6) { if (sigabc.size() < 6) {
sigabc.clear(); sigabc.clear();
return result; return result;
...@@ -611,11 +617,10 @@ double IndexingUtils::Optimize_UB(DblMatrix &UB, ...@@ -611,11 +617,10 @@ double IndexingUtils::Optimize_UB(DblMatrix &UB,
for (int c = 0; c < 3; c++) { for (int c = 0; c < 3; c++) {
UB[r][c] += SMALL; UB[r][c] += SMALL;
GetLatticeParameters(UB, latNew); GetLatticeParameters(UB, latNew);
UB[r][c] -= SMALL; UB[r][c] -= SMALL;
for (int l = 0; l < 7; l++) for (size_t l = 0; l < 7; l++)
derivs[c][l] = (latNew[l] - latOrig[l]) / SMALL; derivs[c][l] = (latNew[l] - latOrig[l]) / SMALL;
} }
...@@ -880,12 +885,8 @@ double IndexingUtils::Optimize_Direction(V3D &best_vec, ...@@ -880,12 +885,8 @@ double IndexingUtils::Optimize_Direction(V3D &best_vec,
@param UB This will be set to the UB matrix that best @param UB This will be set to the UB matrix that best
indexes the supplied list of q_vectors. indexes the supplied list of q_vectors.
@param q_vectors List of locations of peaks in "Q". @param q_vectors List of locations of peaks in "Q".
@param a Lattice parameter "a". @param cell Unit cell defining the parameters a,b,c and
@param b Lattice parameter "b". alpha, beta, gamma.
@param c Lattice parameter "c".
@param alpha Lattice parameter alpha.
@param beta Lattice parameter beta.
@param gamma Lattice parameter gamma.
@param degrees_per_step The number of degrees per step used when @param degrees_per_step The number of degrees per step used when
scanning through all possible directions and scanning through all possible directions and
orientations for the unit cell. NOTE: The orientations for the unit cell. NOTE: The
...@@ -906,14 +907,21 @@ double IndexingUtils::Optimize_Direction(V3D &best_vec, ...@@ -906,14 +907,21 @@ double IndexingUtils::Optimize_Direction(V3D &best_vec,
*/ */
double IndexingUtils::ScanFor_UB(DblMatrix &UB, double IndexingUtils::ScanFor_UB(DblMatrix &UB,
const std::vector<V3D> &q_vectors, double a, const std::vector<V3D> &q_vectors,
double b, double c, double alpha, double beta, const UnitCell &cell, double degrees_per_step,
double gamma, double degrees_per_step,
double required_tolerance) { double required_tolerance) {
if (UB.numRows() != 3 || UB.numCols() != 3) { if (UB.numRows() != 3 || UB.numCols() != 3) {
throw std::invalid_argument("Find_UB(): UB matrix NULL or not 3X3"); throw std::invalid_argument("Find_UB(): UB matrix NULL or not 3X3");
} }
auto a = cell.a();
auto b = cell.b();
auto c = cell.c();
auto alpha = cell.alpha();
auto beta = cell.beta();
auto gamma = cell.gamma();
V3D a_dir; V3D a_dir;
V3D b_dir; V3D b_dir;
V3D c_dir; V3D c_dir;
......
Framework/Geometry/test/IndexingUtilsTest.h
+ 7
− 18
View file @ a5008c56
...@@ -111,12 +111,7 @@ public: ...@@ -111,12 +111,7 @@ public:
std::vector<V3D> q_vectors = getNatroliteQs(); std::vector<V3D> q_vectors = getNatroliteQs();
double a = 6.6; OrientedLattice lattice(6.6, 9.7, 9.9, 84, 71, 70);
double b = 9.7;
double c = 9.9;
double alpha = 84;
double beta = 71;
double gamma = 70;
// test both default case(-1) and // test both default case(-1) and
// case with specified base index(4) // case with specified base index(4)
for (int base_index = -1; base_index < 5; base_index += 5) { for (int base_index = -1; base_index < 5; base_index += 5) {
...@@ -124,9 +119,9 @@ public: ...@@ -124,9 +119,9 @@ public:
size_t num_initial = 3; size_t num_initial = 3;
double degrees_per_step = 3; double degrees_per_step = 3;
double error = IndexingUtils::Find_UB( double error =
UB, q_vectors, a, b, c, alpha, beta, gamma, required_tolerance, IndexingUtils::Find_UB(UB, q_vectors, lattice, required_tolerance,
base_index, num_initial, degrees_per_step); base_index, num_initial, degrees_per_step);
// std::cout << std::endl << "USING LATTICE PARAMETERS\n"; // std::cout << std::endl << "USING LATTICE PARAMETERS\n";
// ShowIndexingStats( UB, q_vectors, required_tolerance ); // ShowIndexingStats( UB, q_vectors, required_tolerance );
...@@ -260,18 +255,12 @@ public: ...@@ -260,18 +255,12 @@ public:
Matrix<double> UB(3, 3, false); Matrix<double> UB(3, 3, false);
int degrees_per_step = 3; int degrees_per_step = 3;
double required_tolerance = 0.2; double required_tolerance = 0.2;
double a = 6.6f;
double b = 9.7f;
double c = 9.9f;
double alpha = 84;
double beta = 71;
double gamma = 70;
UnitCell cell(6.6f, 9.7f, 9.9f, 84, 71, 70);
std::vector<V3D> q_vectors = getNatroliteQs(); std::vector<V3D> q_vectors = getNatroliteQs();
double error = double error = IndexingUtils::ScanFor_UB(
IndexingUtils::ScanFor_UB(UB, q_vectors, a, b, c, alpha, beta, gamma, UB, q_vectors, cell, degrees_per_step, required_tolerance);
degrees_per_step, required_tolerance);
TS_ASSERT_DELTA(error, 0.147397, 1.e-5); TS_ASSERT_DELTA(error, 0.147397, 1.e-5);
......
docs/source/release/v3.9.0/diffraction.rst
+ 2
− 0
View file @ a5008c56
...@@ -7,6 +7,8 @@ Diffraction Changes ...@@ -7,6 +7,8 @@ Diffraction Changes
Crystal Improvements Crystal Improvements
-------------------- --------------------
:ref:`algm-FindUBUsingLatticeParameters` will now return an oriented lattice even when the number of peaks used is very low.
:ref:`algm-FindUBUsingLatticeParameters` has a new option to fix lattice parameters. This will find an orientation, but without optimisation between indexed HKLs and q vectors.
Engineering Diffraction Engineering Diffraction
----------------------- -----------------------
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment