diff --git a/Framework/Crystal/src/FindUBUsingLatticeParameters.cpp b/Framework/Crystal/src/FindUBUsingLatticeParameters.cpp index c1c2f8ea137657d2d3a63e45c8760b0ee22ce9d2..ecc9e59e12b6cfaa61b772e495cb611793d71254 100644 --- a/Framework/Crystal/src/FindUBUsingLatticeParameters.cpp +++ b/Framework/Crystal/src/FindUBUsingLatticeParameters.cpp @@ -80,9 +80,10 @@ void FindUBUsingLatticeParameters::exec() { q_vectors.push_back(peaks[i].getQSampleFrame()); Matrix<double> UB(3, 3, false); - double error = IndexingUtils::Find_UB(UB, q_vectors, a, b, c, alpha, beta, - gamma, tolerance, base_index, - num_initial, degrees_per_step, fixAll); + OrientedLattice lattice(a, b, c, alpha, beta, gamma); + double error = IndexingUtils::Find_UB(UB, q_vectors, lattice, tolerance, + base_index, num_initial, + degrees_per_step, fixAll); g_log.notice() << "Error = " << error << '\n'; g_log.notice() << "UB = " << UB << '\n'; @@ -94,20 +95,6 @@ void FindUBUsingLatticeParameters::exec() { g_log.notice(std::string("UB NOT SAVED.")); } else // tell user how many would be indexed { // 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); - -//// if(!fixAll) { -//// IndexingUtils::Optimize_UB(UB, miller_ind, indexed_qs, sigabc, fixAngles); -//// } - char logInfo[200]; int num_indexed = IndexingUtils::NumberIndexed(UB, q_vectors, tolerance); sprintf(logInfo, @@ -117,27 +104,21 @@ void FindUBUsingLatticeParameters::exec() { num_indexed, n_peaks, tolerance); g_log.notice(std::string(logInfo)); - OrientedLattice o_lattice; - o_lattice.setUB(UB); - o_lattice.setError(sigabc[0], sigabc[1], sigabc[2], sigabc[3], sigabc[4], - sigabc[5]); - - o_lattice.setUB(UB); - 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(); + double calc_a = lattice.a(); + double calc_b = lattice.b(); + double calc_c = lattice.c(); + double calc_alpha = lattice.alpha(); + double calc_beta = lattice.beta(); + double calc_gamma = lattice.gamma(); // 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 " "%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_beta - beta, calc_gamma - gamma); g_log.notice(std::string(logInfo)); - ws->mutableSample().setOrientedLattice(&o_lattice); + ws->mutableSample().setOrientedLattice(&lattice); } } diff --git a/Framework/Geometry/inc/MantidGeometry/Crystal/IndexingUtils.h b/Framework/Geometry/inc/MantidGeometry/Crystal/IndexingUtils.h index 67cecfc0fabdf6b00798c97808e2cd3cb7c2bf0e..c280430d48bdf78885dad17a502811b76e4baf3b 100644 --- a/Framework/Geometry/inc/MantidGeometry/Crystal/IndexingUtils.h +++ b/Framework/Geometry/inc/MantidGeometry/Crystal/IndexingUtils.h @@ -7,6 +7,7 @@ // Includes //---------------------------------------------------------------------- #include "MantidGeometry/DllConfig.h" +#include "MantidGeometry/Crystal/OrientedLattice.h" #include "MantidKernel/V3D.h" #include "MantidKernel/Matrix.h" #include "MantidKernel/Logger.h" @@ -51,11 +52,10 @@ public: /// Find the UB matrix that most nearly indexes the specified qxyz values /// given the lattice parameters static double Find_UB(Kernel::DblMatrix &UB, - const std::vector<Kernel::V3D> &q_vectors, double a, - double b, double c, double alpha, double beta, - double gamma, double required_tolerance, int base_index, - size_t num_initial, double degrees_per_step, - bool fixAll = false); + const std::vector<Kernel::V3D> &q_vectors, + OrientedLattice& lattice, double required_tolerance, + int base_index, size_t num_initial, + double degrees_per_step, bool fixAll = false); /// Find the UB matrix that most nearly indexes the specified qxyz values /// given the range of possible real space unit cell edge lengths. @@ -89,9 +89,9 @@ public: /// Scan rotations to find UB that indexes peaks given lattice parameters static double ScanFor_UB(Kernel::DblMatrix &UB, - const std::vector<Kernel::V3D> &q_vectors, double a, - double b, double c, double alpha, double beta, - double gamma, double degrees_per_step, + const std::vector<Kernel::V3D> &q_vectors, + const UnitCell& lattice, + double degrees_per_step, double required_tolerance); /// Get list of possible directions and lengths for real space unit cell diff --git a/Framework/Geometry/src/Crystal/IndexingUtils.cpp b/Framework/Geometry/src/Crystal/IndexingUtils.cpp index f0e49f6e355502ffa63a200640d0fcbeff4e9dca..e81c586843eecf4f83028dbf2775fef1f8f49a8f 100644 --- a/Framework/Geometry/src/Crystal/IndexingUtils.cpp +++ b/Framework/Geometry/src/Crystal/IndexingUtils.cpp @@ -1,6 +1,5 @@ #include "MantidGeometry/Crystal/IndexingUtils.h" #include "MantidGeometry/Crystal/NiggliCell.h" -#include "MantidGeometry/Crystal/OrientedLattice.h" #include "MantidKernel/Quat.h" #include <algorithm> #include <cmath> @@ -53,12 +52,9 @@ const constexpr double RAD_TO_DEG = 180. / M_PI; @param q_vectors std::vector of V3D objects that contains the list of q_vectors that are to be indexed NOTE: There must be at least 2 q_vectors. - @param a First unit cell edge length in Angstroms. - @param b Second unit cell edge length in Angstroms. - @param c Third unit cell edge length in Angstroms. - @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 lattice The orientated lattice with the lattice + parameters a,b,c and alpha, beta, gamma. The found + UB and errors will be set on this lattice. @param required_tolerance The maximum allowed deviation of Miller indices from integer values for a peak to be indexed. @param base_index The sequence number of the peak that should @@ -87,8 +83,7 @@ const constexpr double RAD_TO_DEG = 180. / M_PI; is <= 0. */ double IndexingUtils::Find_UB(DblMatrix &UB, const std::vector<V3D> &q_vectors, - double a, double b, double c, double alpha, - double beta, double gamma, + OrientedLattice& lattice, double required_tolerance, int base_index, size_t num_initial, double degrees_per_step, bool fixAll) { @@ -157,8 +152,7 @@ double IndexingUtils::Find_UB(DblMatrix &UB, const std::vector<V3D> &q_vectors, for (size_t i = 0; i < num_initial; i++) some_qs.push_back(sorted_qs[i]); - ScanFor_UB(UB, some_qs, a, b, c, alpha, beta, gamma, degrees_per_step, - required_tolerance); + ScanFor_UB(UB, some_qs, lattice, degrees_per_step, required_tolerance); double fit_error = 0; @@ -167,19 +161,11 @@ double IndexingUtils::Find_UB(DblMatrix &UB, const std::vector<V3D> &q_vectors, miller_ind.reserve(q_vectors.size()); indexed_qs.reserve(q_vectors.size()); - // If the user has requested we fix all parameters then return here with the - // found UB and do not optimise it. - if (fixAll) { - GetIndexedPeaks(UB, q_vectors, required_tolerance, miller_ind, indexed_qs, - fit_error); - return fit_error; - } - // now gradually bring in the remaining // peaks and re-optimize the UB to index // them as well size_t count = 0; - while (num_initial < sorted_qs.size()) { + while (!fixAll && num_initial < sorted_qs.size()) { count++; num_initial = std::lround(1.5 * static_cast<double>(num_initial + 3)); // add 3, in case we started with @@ -191,25 +177,24 @@ double IndexingUtils::Find_UB(DblMatrix &UB, const std::vector<V3D> &q_vectors, some_qs.push_back(sorted_qs[i]); try { - std::vector<double> sigabc(7); GetIndexedPeaks(UB, some_qs, required_tolerance, miller_ind, indexed_qs, fit_error); Matrix<double> temp_UB(3, 3, false); - fit_error = Optimize_UB(temp_UB, miller_ind, indexed_qs, sigabc); + fit_error = Optimize_UB(temp_UB, miller_ind, indexed_qs); UB = temp_UB; } catch (...) { // failed to fit using these peaks, so add some more and try again } } - 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 { - std::vector<double> sigabc(7); GetIndexedPeaks(UB, q_vectors, required_tolerance, miller_ind, indexed_qs, fit_error); - Matrix<double> temp_UB(3, 3, false); - fit_error = Optimize_UB(temp_UB, miller_ind, indexed_qs, sigabc); + Matrix<double> temp_UB = UB; + fit_error = Optimize_UB(temp_UB, miller_ind, indexed_qs); UB = temp_UB; } catch (...) { // failed to improve UB using these peaks, so just return the current UB @@ -220,6 +205,10 @@ double IndexingUtils::Find_UB(DblMatrix &UB, const std::vector<V3D> &q_vectors, // HKL space. GetIndexedPeaks(UB, q_vectors, required_tolerance, miller_ind, indexed_qs, 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; } @@ -894,12 +883,8 @@ double IndexingUtils::Optimize_Direction(V3D &best_vec, @param UB This will be set to the UB matrix that best indexes the supplied list of q_vectors. @param q_vectors List of locations of peaks in "Q". - @param a Lattice parameter "a". - @param b Lattice parameter "b". - @param c Lattice parameter "c". - @param alpha Lattice parameter alpha. - @param beta Lattice parameter beta. - @param gamma Lattice parameter gamma. + @param cell Unit cell defining the parameters a,b,c and + alpha, beta, gamma. @param degrees_per_step The number of degrees per step used when scanning through all possible directions and orientations for the unit cell. NOTE: The @@ -920,14 +905,22 @@ double IndexingUtils::Optimize_Direction(V3D &best_vec, */ double IndexingUtils::ScanFor_UB(DblMatrix &UB, - const std::vector<V3D> &q_vectors, double a, - double b, double c, double alpha, double beta, - double gamma, double degrees_per_step, + const std::vector<V3D> &q_vectors, + const UnitCell& cell, + double degrees_per_step, double required_tolerance) { if (UB.numRows() != 3 || UB.numCols() != 3) { 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 b_dir; V3D c_dir; diff --git a/Framework/Geometry/test/IndexingUtilsTest.h b/Framework/Geometry/test/IndexingUtilsTest.h index b8b66458b2b4b6dfe85472af9af4e42ff788d271..0974588a716466556c90e4c07104ddbaaa7e5740 100644 --- a/Framework/Geometry/test/IndexingUtilsTest.h +++ b/Framework/Geometry/test/IndexingUtilsTest.h @@ -111,12 +111,7 @@ public: std::vector<V3D> q_vectors = getNatroliteQs(); - double a = 6.6; - double b = 9.7; - double c = 9.9; - double alpha = 84; - double beta = 71; - double gamma = 70; + OrientedLattice lattice(6.6, 9.7, 9.9, 84, 71, 70); // test both default case(-1) and // case with specified base index(4) for (int base_index = -1; base_index < 5; base_index += 5) { @@ -125,7 +120,7 @@ public: double degrees_per_step = 3; double error = IndexingUtils::Find_UB( - UB, q_vectors, a, b, c, alpha, beta, gamma, required_tolerance, + UB, q_vectors, lattice, required_tolerance, base_index, num_initial, degrees_per_step); // std::cout << std::endl << "USING LATTICE PARAMETERS\n"; @@ -260,17 +255,12 @@ public: Matrix<double> UB(3, 3, false); int degrees_per_step = 3; 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(); double error = - IndexingUtils::ScanFor_UB(UB, q_vectors, a, b, c, alpha, beta, gamma, + IndexingUtils::ScanFor_UB(UB, q_vectors, cell, degrees_per_step, required_tolerance); TS_ASSERT_DELTA(error, 0.147397, 1.e-5);