Skip to content

GitLab

Commit 6a23a9f1 authored by Martyn Gigg's avatar Martyn Gigg
Browse files

Significant refactor of IndexPeaks 

Improves readability, code reuse and style.

Refs #26829
parent eac26e44
Expand all Hide whitespace changes
Inline Side-by-side
Framework/Crystal/inc/MantidCrystal/IndexPeaks.h
View file @ 6a23a9f1
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2011 ISIS Rutherford Appleton Laboratory UKRI,
// Copyright © 2019 ISIS Rutherford Appleton Laboratory UKRI,
// NScD Oak Ridge National Laboratory, European Spallation Source
// & Institut Laue - Langevin
// SPDX - License - Identifier: GPL - 3.0 +
#ifndef MANTID_CRYSTAL_INDEX_PEAKS_H_
#define MANTID_CRYSTAL_INDEX_PEAKS_H_
#include "MantidAPI/Algorithm.h"
#include "MantidKernel/System.h"
namespace Mantid {
namespace Crystal {
/** IndexPeaks : Algorithm to use the UB saved in the sample associated
with the specified PeaksWorkspace, to index the peaks in the workspace.
@author Dennis Mikkelson
@date 2011-08-17
*/
/**
* Implements an algorithm for indexing main and satellites peaks
* in single crystal peaks.
*/
class DLLExport IndexPeaks : public API::Algorithm {
public:
/// Algorithm's name for identification
const std::string name() const override { return "IndexPeaks"; };
/// Summary of algorithms purpose
const std::string name() const override { return "IndexPeaks"; }
const std::string summary() const override {
return "Index the peaks using the UB from the sample.";
}
/// Algorithm's version for identification
int version() const override { return 1; }
const std::vector<std::string> seeAlso() const override {
return {"IndexSXPeaks"};
}
/// Algorithm's category for identification
const std::string category() const override { return "Crystal\\Peaks"; }
std::map<std::string, std::string> validateInputs() override;
private:
/// Initialise the properties
protected:
void init() override;
/// Run the algorithm
void exec() override;
};
......
Framework/Crystal/src/IndexPeaks.cpp
View file @ 6a23a9f1
This diff is collapsed.
Framework/Crystal/test/IndexPeaksTest.h
View file @ 6a23a9f1
......@@ -130,33 +130,62 @@ public:
assertNumberPeaksIndexed(*alg, 0, 0, 0);
}
void test_no_roundHKLS_leaves_peaks_as_found() {
void test_no_commonUB_optimizes_UB_per_run_for_no_satellites() {
const auto ws = createTestPeaksWorkspaceMainReflOnly();
const auto mainToleranceAsStr = "0.1";
auto alg =
indexPeaks(ws, {{"Tolerance", mainToleranceAsStr}, {"RoundHKLs", "0"}});
// Check that the peaks were all indexed
const double mainTolerance{std::stod(mainToleranceAsStr)};
auto &peaks = ws->getPeaks();
for (const auto &peak : peaks) {
TS_ASSERT(IndexingUtils::ValidIndex(peak.getHKL(), mainTolerance))
}
// Change some run numbers
ws->getPeak(0).setRunNumber(3008);
ws->getPeak(4).setRunNumber(3008);
auto alg = indexPeaks(
ws,
{{"Tolerance", "0.1"}, {"RoundHKLs", "1"}, {"CommonUBForAll", "0"}});
// Check the output properties
assertNumberPeaksIndexed(*alg, 5, 5, 0);
const double averageError = alg->getProperty("AverageError");
TS_ASSERT_DELTA(averageError, 0.00639, 1e-4)
assertErrorsAsExpected(*alg, 0.00883, 0.00883, 0.);
// spot check a few peaks for
// fractional Miller indices
V3D peak_0_hkl(4, 1, 6);
V3D peak_1_hkl(3, -1, 4);
V3D peak_2_hkl(4, -1, 5);
V3D peak_3_hkl(3, -0, 7);
V3D peak_4_hkl(2, -4, 7);
const auto &peaks = ws->getPeaks();
V3D error = peak_0_hkl - peaks[0].getHKL();
TS_ASSERT_DELTA(error.norm(), 0.0, 2e-4)
error = peak_1_hkl - peaks[1].getHKL();
TS_ASSERT_DELTA(error.norm(), 0.0, 1e-4)
error = peak_2_hkl - peaks[2].getHKL();
TS_ASSERT_DELTA(error.norm(), 0.0, 1e-4)
error = peak_3_hkl - peaks[3].getHKL();
TS_ASSERT_DELTA(error.norm(), 0.0, 1e-4)
error = peak_4_hkl - peaks[4].getHKL();
TS_ASSERT_DELTA(error.norm(), 0.0, 2e-4)
}
void test_no_roundHKLS_leaves_peaks_as_found() {
const auto ws = createTestPeaksWorkspaceMainReflOnly();
auto alg = indexPeaks(ws, {{"Tolerance", "0.1"}, {"RoundHKLs", "0"}});
// Check the output properties
assertNumberPeaksIndexed(*alg, 5, 5, 0);
assertErrorsAsExpected(*alg, 0.00639, 0.00639, 0.);
// spot check a few peaks for
// fractional Miller indices
const V3D peak_0_hkl_d(4.00682, 0.97956, 5.99368); // first peak
const V3D peak_1_hkl_d(2.99838, -0.99760, 4.00141);
const V3D peak_2_hkl_d(3.99737, -0.99031, 5.00250);
const V3D peak_3_hkl_d(2.99419, 0.01736, 7.00538);
const V3D peak_4_hkl_d(2.00277, -4.00813, 6.99744); // last peak
const auto &peaks = ws->getPeaks();
V3D error = peak_0_hkl_d - peaks[0].getHKL();
TS_ASSERT_DELTA(error.norm(), 0.0, 2e-4)
......@@ -185,16 +214,9 @@ public:
TS_ASSERT(IndexingUtils::ValidIndex(peak.getHKL(), mainTolerance))
}
// Check that the peaks were all indexed
for (const auto &peak : peaks) {
TS_ASSERT(IndexingUtils::ValidIndex(peak.getHKL(), mainTolerance))
}
// Check the output properties
assertNumberPeaksIndexed(*alg, 5, 5, 0);
const double averageError = alg->getProperty("AverageError");
TS_ASSERT_DELTA(averageError, 0.00673, 1e-3)
assertErrorsAsExpected(*alg, 0.00639, 0.00639, 0.);
// spot check a few peaks for
// integer Miller indices
......@@ -221,7 +243,7 @@ public:
}
void
test_zero_satellite_tol_only_indexes_main_refl_with_modvectors_on_from_ub() {
test_zero_satellite_tol_only_indexes_main_refl_with_modvectors_from_ub() {
const auto peaksWS =
createTestPeaksWorkspaceWithSatellites(1, {V3D(0.333, -0.667, 0.333)});
const auto mainTolerance{"0.1"}, sateTolerance{"0."};
......@@ -231,6 +253,10 @@ public:
{"ToleranceForSatellite", sateTolerance}});
assertNumberPeaksIndexed(*alg, 2, 2, 0);
std::cerr << "\n";
for (int i = 0; i < peaksWS->getNumberPeaks(); ++i) {
std::cerr << peaksWS->getPeak(i).getHKL() << "\n";
}
assertIndexesAsExpected(*peaksWS,
{V3D(-1, 2, -9), V3D(0, 0, 0), V3D(-1, 1, -10),
V3D(0, 0, 0), V3D(0, 0, 0)});
......@@ -251,7 +277,7 @@ public:
V3D(0, 0, 0), V3D(0, 0, 0)});
}
void test_exec_with_common_ub_indexes_main_reflections_only() {
void test_exec_with_common_ub_and_zero_mod_vectors_indexes_main_refl_only() {
const auto peaksWS =
createTestPeaksWorkspaceWithSatellites(0, std::vector<V3D>());
const auto sateTolerance{"1."};
......@@ -267,6 +293,22 @@ public:
assertErrorsAsExpected(*alg, 0.0140447, 0.0140447, 0.);
}
void xtest_exec_with_common_ub_and_non_zero_mod_vectors_indexes_both() {
const auto peaksWS =
createTestPeaksWorkspaceWithSatellites(1, {V3D(0.333, -0.667, 0.333)});
const auto sateTolerance{"1."};
const auto alg =
indexPeaks(peaksWS, {{"CommonUBForAll", "1"},
{"ToleranceForSatellite", sateTolerance}});
assertNumberPeaksIndexed(*alg, 5, 2, 3);
assertIndexesAsExpected(*peaksWS,
{V3D(-1, 2, -9), V3D(0, 0, 0), V3D(-1, 1, -10),
V3D(0, 0, 0), V3D(0, 0, 0)});
assertErrorsAsExpected(*alg, 0.0140447, 0.0140447, 0.);
}
// --------------------------- Failure tests -----------------------------
void test_workspace_with_no_oriented_lattice_gives_validateInputs_error() {
......
Framework/DataObjects/inc/MantidDataObjects/Peak.h
View file @ 6a23a9f1
......@@ -104,6 +104,7 @@ public:
double getK() const override;
double getL() const override;
Mantid::Kernel::V3D getHKL() const override;
bool isIndexed() const override;
Mantid::Kernel::V3D getIntHKL() const override;
Mantid::Kernel::V3D getIntMNP() const override;
void setH(double m_H) override;
......
Framework/DataObjects/src/Peak.cpp
View file @ 6a23a9f1
......@@ -759,6 +759,13 @@ double Peak::getL() const { return m_L; }
/** Return the HKL vector */
Mantid::Kernel::V3D Peak::getHKL() const { return V3D(m_H, m_K, m_L); }
/** Return True if the peak has been indexed */
bool Peak::isIndexed() const {
if (m_H == 0. && m_K == 0. && m_L == 0.)
return false;
return true;
}
/** Return the int HKL vector */
Mantid::Kernel::V3D Peak::getIntHKL() const { return m_intHKL; }
......
Framework/DataObjects/test/PeakTest.h
View file @ 6a23a9f1
......@@ -236,6 +236,13 @@ public:
TS_ASSERT_EQUALS(p.getHKL(), V3D(1.0, 2.0, 3.0));
}
void test_isIndexed() {
Peak p(inst, 10000, 2.0);
TS_ASSERT_EQUALS(false, p.isIndexed());
p.setHKL(1, 2, 3);
TS_ASSERT_EQUALS(true, p.isIndexed());
}
void test_samplePos() {
Peak p(inst, 10000, 2.0);
p.setSamplePos(1.0, 1.0, 1.0);
......
Framework/Geometry/inc/MantidGeometry/Crystal/IPeak.h
View file @ 6a23a9f1
......@@ -45,6 +45,7 @@ public:
virtual double getK() const = 0;
virtual double getL() const = 0;
virtual Mantid::Kernel::V3D getHKL() const = 0;
virtual bool isIndexed() const = 0;
virtual Mantid::Kernel::V3D getIntHKL() const = 0;
virtual void setH(double m_H) = 0;
virtual void setK(double m_K) = 0;
......
Framework/Geometry/inc/MantidGeometry/Crystal/IndexingUtils.h
View file @ 6a23a9f1
......@@ -194,6 +194,12 @@ public:
std::vector<Kernel::V3D> &miller_indices,
double &ave_error);
/// Given a UB, calculate the miller indices for given q vector
static bool CalculateMillerIndices(const Kernel::DblMatrix &inverseUB,
const Kernel::V3D &q_vector,
double tolerance,
Kernel::V3D &miller_indices);
/// Get lists of indices and Qs for peaks indexed in the specified direction
static int GetIndexedPeaks_1D(const Kernel::V3D &direction,
const std::vector<Kernel::V3D> &q_vectors,
......
Framework/Geometry/src/Crystal/IndexingUtils.cpp
View file @ 6a23a9f1
......@@ -2420,20 +2420,14 @@ int IndexingUtils::CalculateMillerIndices(const DblMatrix &UB,
miller_indices.reserve(q_vectors.size());
int count = 0;
double h_error, k_error, l_error;
ave_error = 0.0;
V3D hkl;
for (const auto &q_vector : q_vectors) {
hkl = UB_inverse * q_vector / (2.0 * M_PI);
if (ValidIndex(hkl, tolerance)) {
V3D hkl;
if (CalculateMillerIndices(UB_inverse, q_vector, tolerance, hkl)) {
count++;
miller_indices.emplace_back(hkl);
h_error = std::abs(std::round(hkl[0]) - hkl[0]);
k_error = std::abs(std::round(hkl[1]) - hkl[1]);
l_error = std::abs(std::round(hkl[2]) - hkl[2]);
ave_error += h_error + k_error + l_error;
} else
miller_indices.emplace_back(0, 0, 0);
ave_error += hkl.hklError();
}
miller_indices.emplace_back(std::move(hkl));
}
if (count > 0) {
......@@ -2443,6 +2437,38 @@ int IndexingUtils::CalculateMillerIndices(const DblMatrix &UB,
return count;
}
/**
Calculate the Miller Indices for each of the specified Q vectors, using the
inverse of the specified UB matrix. If the peaks could not be indexed it is
set to (0,0,0)
@param UB A 3x3 matrix of doubles holding the inverse UB matrix.
The matrix is not checked for validity
@param q_vectors std::vector of V3D objects that contains the list of
q_vectors that are to be indexed.
@param tolerance The maximum allowed distance between each component
of UB^(-1)*Q and the nearest integer value, required to
to count the peak as indexed by UB.
@param miller_indices This vector returns a list of Miller Indices, with
one entry for each given Q vector.
@return True if the peak was index, false otherwise
*/
bool IndexingUtils::CalculateMillerIndices(const DblMatrix &inverseUB,
const V3D &q_vector,
double tolerance,
V3D &miller_indices) {
miller_indices = inverseUB * q_vector / (2.0 * M_PI);
if (ValidIndex(miller_indices, tolerance)) {
return true;
} else {
miller_indices = V3D(0, 0, 0);
return false;
}
}
/**
Given one plane normal direction for a family of parallel planes in
reciprocal space, find the peaks that lie on these planes to within the
......@@ -2515,7 +2541,8 @@ int IndexingUtils::GetIndexedPeaks_1D(const V3D &direction,
to the unit cell edges, then the resulting indices will be proper Miller
indices for the peaks. This method is similar to GetIndexedPeaks_1D, but
checks three directions simultaneously and requires that the peak lies
on all three families of planes simultaneously and does NOT index as (0,0,0).
on all three families of planes simultaneously and does NOT index as
(0,0,0).
@param direction_1 Direction vector in the direction of the normal
vector for the first family of parallel planes.
......@@ -2786,16 +2813,15 @@ std::vector<V3D> IndexingUtils::MakeCircleDirections(int n_steps,
plane_spacing. The direction is chosen from the specified direction_list.
@param best_direction This will be set to the direction that minimizes
the sum squared distances of projections of peaks
from integer multiples of the specified plane
spacing.
the sum squared distances of projections of
peaks from integer multiples of the specified plane spacing.
@param q_vectors List of peak positions, specified according to
the convention that |q| = 1/d. (i.e. Q/2PI)
@param direction_list List of possible directions for plane normals.
Initially, this will be a long list of possible
directions from MakeHemisphereDirections().
@param plane_spacing The required spacing between planes in reciprocal
space.
@param plane_spacing The required spacing between planes in
reciprocal space.
@param required_tolerance The maximum deviation of the component of a
peak Qxyz in the direction of the best_direction
vector for that peak to count as being indexed.
......@@ -2863,7 +2889,8 @@ int IndexingUtils::SelectDirection(V3D &best_direction,
* @param UB A non-singular matrix representing an orientation
* matrix.
* @param lattice_par std::vector of doubles that will contain the lattice
* parameters and cell volume as it's first seven entries.
* parameters and cell volume as it's first seven
* entries.
* @return true if the lattice_par vector was filled with the lattice
* parameters and false if the matrix could not be inverted.
*/
......
Framework/Geometry/test/IndexingUtilsTest.h
View file @ 6a23a9f1
......@@ -637,6 +637,22 @@ public:
}
}
void test_CalculateMillerIndicesSingleQ_for_Valid_Index() {
const auto q_vectors = getNatroliteQs();
const auto indices = getNatroliteIndices();
auto UB = getNatroliteUB();
UB.Invert();
const double tolerance = 0.1;
V3D miller_indices;
const bool success = IndexingUtils::CalculateMillerIndices(
UB, q_vectors[0], tolerance, miller_indices);
TS_ASSERT(success);
const auto diff = (indices[0] - miller_indices).norm();
TS_ASSERT_DELTA(diff, 0, 0.1);
}
void test_GetIndexedPeaks_1D() {
int correct_indices[] = {1, 4, 2, 0, 1, 3, 0, -1, 0, -1, -2, -3};
......
Framework/Geometry/test/MockObjects.h
View file @ 6a23a9f1
......@@ -83,6 +83,7 @@ public:
MOCK_CONST_METHOD0(getK, double());
MOCK_CONST_METHOD0(getL, double());
MOCK_CONST_METHOD0(getHKL, Mantid::Kernel::V3D());
MOCK_CONST_METHOD0(isIndexed, bool());
MOCK_CONST_METHOD0(getIntHKL, Mantid::Kernel::V3D());
MOCK_CONST_METHOD0(getSamplePos, Mantid::Kernel::V3D());
MOCK_METHOD1(setH, void(double m_H));
......@@ -90,7 +91,7 @@ public:
MOCK_METHOD1(setL, void(double m_L));
MOCK_METHOD3(setHKL, void(double H, double K, double L));
MOCK_METHOD1(setHKL, void(const Mantid::Kernel::V3D &HKL));
MOCK_METHOD1(setIntHKL, void(const Mantid::Kernel::V3D HKL));
MOCK_METHOD1(setIntHKL, void(const Mantid::Kernel::V3D &HKL));
MOCK_METHOD3(setSamplePos, void(double samX, double samY, double samZ));
MOCK_METHOD1(setSamplePos, void(const Mantid::Kernel::V3D &XYZ));
MOCK_CONST_METHOD0(getQLabFrame, Mantid::Kernel::V3D());
......
Markdown is supported
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