From 6a23a9f16ab5e6196d6c3cc4738138bf3ec2508d Mon Sep 17 00:00:00 2001
From: Martyn Gigg <martyn.gigg@gmail.com>
Date: Tue, 17 Sep 2019 16:46:14 +0100
Subject: [PATCH] Significant refactor of IndexPeaks
Improves readability, code reuse and style.
Refs #26829
---
.../Crystal/inc/MantidCrystal/IndexPeaks.h | 27 +-
Framework/Crystal/src/IndexPeaks.cpp | 1143 +++++++++++------
Framework/Crystal/test/IndexPeaksTest.h | 88 +-
.../DataObjects/inc/MantidDataObjects/Peak.h | 1 +
Framework/DataObjects/src/Peak.cpp | 7 +
Framework/DataObjects/test/PeakTest.h | 7 +
.../inc/MantidGeometry/Crystal/IPeak.h | 1 +
.../MantidGeometry/Crystal/IndexingUtils.h | 6 +
.../Geometry/src/Crystal/IndexingUtils.cpp | 63 +-
Framework/Geometry/test/IndexingUtilsTest.h | 16 +
Framework/Geometry/test/MockObjects.h | 3 +-
11 files changed, 906 insertions(+), 456 deletions(-)
diff --git a/Framework/Crystal/inc/MantidCrystal/IndexPeaks.h b/Framework/Crystal/inc/MantidCrystal/IndexPeaks.h
index 897eb610a42..9229b551032 100644
--- a/Framework/Crystal/inc/MantidCrystal/IndexPeaks.h
+++ b/Framework/Crystal/inc/MantidCrystal/IndexPeaks.h
@@ -1,46 +1,35 @@
// 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;
};
diff --git a/Framework/Crystal/src/IndexPeaks.cpp b/Framework/Crystal/src/IndexPeaks.cpp
index 6cf2e353ab9..ea1e7695f4e 100644
--- a/Framework/Crystal/src/IndexPeaks.cpp
+++ b/Framework/Crystal/src/IndexPeaks.cpp
@@ -11,442 +11,795 @@
#include "MantidGeometry/Crystal/OrientedLattice.h"
#include "MantidKernel/BoundedValidator.h"
+#include <boost/optional/optional.hpp>
+
+#include <algorithm>
+
namespace Mantid {
namespace Crystal {
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(IndexPeaks)
-using namespace Mantid::Kernel;
-using namespace Mantid::API;
-using namespace Mantid::DataObjects;
-using namespace Mantid::Geometry;
-
-/** Initialize the algorithm's properties.
- */
-void IndexPeaks::init() {
- this->declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace>>(
- "PeaksWorkspace", "", Direction::InOut),
- "Input Peaks Workspace");
-
- auto mustBePositive = boost::make_shared<BoundedValidator<double>>();
- mustBePositive->setLower(0.0);
-
- this->declareProperty(
- std::make_unique<PropertyWithValue<double>>(
- "Tolerance", 0.15, mustBePositive, Direction::Input),
- "Indexing Tolerance (0.15)");
-
- this->declareProperty(
- std::make_unique<PropertyWithValue<double>>(
- "ToleranceForSatellite", 0.15, mustBePositive, Direction::Input),
- "Satellite Indexing Tolerance (0.15)");
-
- this->declareProperty("RoundHKLs", true,
- "Round H, K and L values to integers");
-
- this->declareProperty("CommonUBForAll", false,
- "Index all orientations with a common UB");
+using DataObjects::Peak;
+using DataObjects::PeaksWorkspace;
+using DataObjects::PeaksWorkspace_sptr;
+using Geometry::IndexingUtils;
+using Geometry::OrientedLattice;
+using Kernel::DblMatrix;
+using Kernel::V3D;
+
+namespace {
+const auto OPTIMIZE_UB_ATTEMPTS{4};
+
+namespace Prop {
+const std::string PEAKSWORKSPACE{"PeaksWorkspace"};
+const std::string TOLERANCE{"Tolerance"};
+const std::string SATE_TOLERANCE{"ToleranceForSatellite"};
+const std::string ROUNDHKLS{"RoundHKLs"};
+const std::string COMMONUB{"CommonUBForAll"};
+const std::string AVERAGE_ERR{"AverageError"};
+const std::string NUM_INDEXED{"NumIndexed"};
+const std::string MAIN_NUM_INDEXED{"MainNumIndexed"};
+const std::string SATE_NUM_INDEXED{"SateNumIndexed"};
+const std::string MAIN_ERR{"MainError"};
+const std::string SATE_ERR{"SatelliteError"};
+
+struct SatelliteIndexingArgs {
+ const double tolerance;
+ const int maxOrder;
+ const std::vector<V3D> modVectors;
+};
+
+struct IndexPeaksArgs {
+ static IndexPeaksArgs parse(const API::Algorithm &alg) {
+ PeaksWorkspace_sptr peaksWS = alg.getProperty(PEAKSWORKSPACE);
+ const auto &lattice = peaksWS->sample().getOrientedLattice();
+
+ auto addIfNonZero = [](const V3D &modVec, std::vector<V3D> &modVectors) {
+ if (std::fabs(modVec[0]) > 0 && std::fabs(modVec[1]) > 0 &&
+ std::fabs(modVec[2]) > 0)
+ modVectors.emplace_back(modVec);
+ };
+ const int maxOrder = lattice.getMaxOrder();
+ std::vector<V3D> modVectorsToUse;
+ if (maxOrder > 0) {
+ modVectorsToUse.reserve(3);
+ for (auto i = 0; i < 3; ++i) {
+ addIfNonZero(lattice.getModVec(i), modVectorsToUse);
+ }
+ }
- this->declareProperty(std::make_unique<PropertyWithValue<double>>(
- "AverageError", 0.0, Direction::Output),
- "Gets set with the average HKL indexing error.");
+ return {peaksWS, alg.getProperty(TOLERANCE), alg.getProperty(ROUNDHKLS),
+ alg.getProperty(COMMONUB),
+ SatelliteIndexingArgs{alg.getProperty(SATE_TOLERANCE), maxOrder,
+ modVectorsToUse}};
+ }
- this->declareProperty(std::make_unique<PropertyWithValue<int>>(
- "NumIndexed", 0, Direction::Output),
- "Gets set with the number of indexed peaks.");
+ PeaksWorkspace_sptr workspace;
+ const double mainTolerance;
+ const bool roundHKLs;
+ const bool commonUB;
+ const SatelliteIndexingArgs satellites;
+};
+} // namespace Prop
- this->declareProperty(std::make_unique<PropertyWithValue<int>>(
- "MainNumIndexed", 0, Direction::Output),
- "Gets set with the number of indexed main peaks.");
+/**
+ * Track details about the peaks successfully indexed
+ */
+struct PeakIndexingStats {
+ PeakIndexingStats &operator+=(const PeakIndexingStats &rhs) {
+ numIndexed += rhs.numIndexed;
+ error += rhs.error;
+ return *this;
+ }
+ int numIndexed{0};
+ double error{0.0};
+};
- this->declareProperty(std::make_unique<PropertyWithValue<int>>(
- "SateNumIndexed", 0, Direction::Output),
- "Gets set with the number of indexed main peaks.");
+/**
+ * Track details of main and satellite reflections that have
+ * been indexed
+ */
+struct CombinedIndexingStats {
+ CombinedIndexingStats &operator+=(const CombinedIndexingStats &rhs) {
+ main += rhs.main;
+ satellites += rhs.satellites;
+ return *this;
+ }
- this->declareProperty(
- std::make_unique<PropertyWithValue<double>>("MainError", 0.0,
- Direction::Output),
- "Gets set with the average HKL indexing error of Main Peaks.");
+ /// Return the total number of peaks indexed
+ inline int totalNumIndexed() const {
+ return main.numIndexed + satellites.numIndexed;
+ }
+ /// Return the number of fundamental peaks indexed
+ inline double mainError() const {
+ if (main.numIndexed == 0)
+ return 0.0;
+ return main.error / main.numIndexed;
+ }
+ /// Return the number of satellite peaks indexed
+ inline double satelliteError() const {
+ if (satellites.numIndexed == 0)
+ return 0.0;
+ return satellites.error / satellites.numIndexed;
+ }
+ /// Return the average error for both main/satellites indexed
+ inline double averageError() const {
+ return (main.error + satellites.error) / totalNumIndexed();
+ }
- this->declareProperty(
- std::make_unique<PropertyWithValue<double>>("SatelliteError", 0.0,
- Direction::Output),
- "Gets set with the average HKL indexing error of Satellite Peaks.");
-}
+ PeakIndexingStats main;
+ PeakIndexingStats satellites;
+};
-/** Execute the algorithm.
+/**
+ * Attempt to optimize the UB for the given set of peaks
+ * @param ubOrig Original sample UB matrix
+ * @param qSample The Q_sample of each peak to optimize
+ * @return A new optimized UB
*/
-void IndexPeaks::exec() {
- PeaksWorkspace_sptr ws = this->getProperty("PeaksWorkspace");
- OrientedLattice o_lattice = ws->mutableSample().getOrientedLattice();
- const Matrix<double> &UB = o_lattice.getUB();
-
- if (!IndexingUtils::CheckUB(UB)) {
- throw std::runtime_error(
- "ERROR: The stored UB is not a valid orientation matrix");
+DblMatrix optimizeUBMatrix(const DblMatrix &ubOrig,
+ const std::vector<V3D> &qSample,
+ const double tolerance) {
+ DblMatrix optimizedUB(ubOrig);
+
+ double errorAtStart{0.0};
+ std::vector<V3D> millerIndices;
+ millerIndices.reserve(qSample.size());
+ const int numIndexedAtStart = IndexingUtils::CalculateMillerIndices(
+ optimizedUB, qSample, tolerance, millerIndices, errorAtStart);
+
+ if (numIndexedAtStart < 3) {
+ // can't optimize without at least 3 indexed peaks
+ return optimizedUB;
}
- bool round_hkls = this->getProperty("RoundHKLs");
- bool commonUB = this->getProperty("CommonUBForAll");
-
- std::vector<Peak> &peaks = ws->getPeaks();
- size_t n_peaks = ws->getNumberPeaks();
- int total_indexed = 0;
- int total_main = 0;
- int total_sate = 0;
- double average_error;
- double average_main_error;
- double average_sate_error;
- double tolerance = this->getProperty("Tolerance");
-
- if (commonUB && o_lattice.getModVec(0) == V3D(0, 0, 0)) {
- std::vector<V3D> miller_indices;
- std::vector<V3D> q_vectors;
-
- q_vectors.reserve(n_peaks);
- for (const auto &peak : peaks) {
- q_vectors.push_back(peak.getQSampleFrame());
+ for (auto i = 0; i < OPTIMIZE_UB_ATTEMPTS; ++i) {
+ try {
+ // optimization requires rounded indices
+ IndexingUtils::RoundHKLs(millerIndices);
+ IndexingUtils::Optimize_UB(optimizedUB, millerIndices, qSample);
+ } catch (...) {
+ // If there is any problem, such as too few
+ // independent peaks, just use the original UB
+ optimizedUB = ubOrig;
+ break;
}
+ double errorInLoop{0.0};
+ const int numIndexedInLoop = IndexingUtils::CalculateMillerIndices(
+ optimizedUB, qSample, tolerance, millerIndices, errorInLoop);
+ if (numIndexedInLoop < numIndexedAtStart) // use the original UB
+ break;
+ }
+ return optimizedUB;
+}
- total_indexed = IndexingUtils::CalculateMillerIndices(
- UB, q_vectors, tolerance, miller_indices, average_error);
- if (round_hkls)
- IndexingUtils::RoundHKLs(miller_indices);
-
- for (size_t i = 0; i < n_peaks; i++) {
- peaks[i].setHKL(miller_indices[i]);
- peaks[i].setIntHKL(miller_indices[i]);
- peaks[i].setIntMNP(V3D(0, 0, 0));
- }
- } else {
- double total_error = 0;
- double total_main_error = 0;
- double total_sate_error = 0;
- double satetolerance = this->getProperty("ToleranceForSatellite");
-
- // get list of run numbers in this peaks workspace
- std::vector<int> run_numbers;
- for (const auto &peak : peaks) {
- int run = peak.getRunNumber();
- bool found = false;
- size_t k = 0;
- while (k < run_numbers.size() && !found) {
- if (run == run_numbers[k])
- found = true;
- else
- k++;
+/// <IntHKL, IntMNP, error>
+using IndexedSatelliteInfo = std::tuple<V3D, V3D, double>;
+
+/**
+ * @brief Attempt to index a satellite reflection given a HKL from a failed
+ * indexing of a main reflection.
+ * - loop over [-maxOrder, maxOrder] (skipping 0) and for each order:
+ * - compute "hkl_main - order x mod_vector" and see if
+ * this is a valid set of indices
+ * - if it is accept it as a satellite
+ * - if not do nothing
+ * this continues for each modulation provided and the satellite will be the
+ * last found with a valid index.
+
+ * @param mainHKL The HKL of a failed attempt to index a fundamental
+ reflection
+ * @param maxOrder The highest order of to search away from the peak
+ * @param modVectors
+ * @param tolerance
+ * @return
+ */
+boost::optional<IndexedSatelliteInfo>
+indexSatellite(const V3D &mainHKL, int maxOrder,
+ const std::vector<V3D> &modVectors, double tolerance) {
+ assert(modVectors.size() <= 3);
+
+ bool foundSatellite{false};
+ V3D candidateIntHKL, candidateMNP;
+ for (auto i = 0u; i < modVectors.size(); ++i) {
+ const auto &modVector = modVectors[i];
+ for (int order = -maxOrder; order <= maxOrder; ++order) {
+ if (order == 0)
+ continue;
+ candidateIntHKL = mainHKL - modVector * order;
+ if (IndexingUtils::ValidIndex(candidateIntHKL, tolerance)) {
+ candidateMNP[i] = order;
+ foundSatellite = true;
+ // we deliberately don't break and use the last valid
+ // reflection we find.
}
- if (!found)
- run_numbers.push_back(run);
}
+ }
+ if (foundSatellite)
+ return std::make_tuple(candidateIntHKL, candidateMNP,
+ candidateIntHKL.hklError());
+ else
+ return boost::none;
+}
- // index the peaks for each run separately, using a UB matrix optimized for
- // that run
+/**
+ * Attempt to index any unindexed peaks as if they were satellites and set
+ * the HKLs of any indexed peaks.
+ *
+ * @param peaks A list of peaks, with some possibly indexed as main
+ * reflections
+ * @param millerIndices A list of miller indices calculated from attempting
+ * to index the peaks in the peaks list as fundamental peaks. The size is
+ * assumed to match the size of peaks.
+ * @param roundHKLs If true round the output HKL indices
+ * @param args Additional arguments relating to the satellites peaks
+ * @returns A PeakIndexingStats object describing the output
+ */
+PeakIndexingStats indexSatellites(const std::vector<Peak *> &peaks,
+ std::vector<V3D> &millerIndices,
+ const bool roundHKLs,
+ const Prop::SatelliteIndexingArgs &args) {
+ PeakIndexingStats stats;
+
+ for (auto i = 0u; i < peaks.size(); ++i) {
+ auto peak = peaks[i];
+ auto &mainHKL{millerIndices[i]};
+
+ // TODO: THIS TEST IS WRONG AS WE HAVE NOT SET THE HKL YET!!!!
+ if (peak->isIndexed()) {
+ if (roundHKLs)
+ IndexingUtils::RoundHKL(mainHKL);
+ peak->setHKL(mainHKL);
+ peak->setIntHKL(mainHKL);
+ peak->setIntMNP(V3D(0, 0, 0));
+ } else {
+ auto result = indexSatellite(mainHKL, args.maxOrder, args.modVectors,
+ args.tolerance);
+ if (result) {
+ const auto &satelliteInfo = result.get();
+ if (roundHKLs)
+ IndexingUtils::RoundHKL(mainHKL);
+ peak->setHKL(mainHKL);
+ peak->setIntHKL(std::get<0>(satelliteInfo));
+ peak->setIntMNP(std::get<1>(satelliteInfo));
+ stats.numIndexed++;
+ stats.error += std::get<2>(satelliteInfo) / 3.;
+ }
+ }
+ }
+ return stats;
+}
- for (const int run : run_numbers) {
- std::vector<V3D> miller_indices;
- std::vector<V3D> q_vectors;
+/**
+ * Index the main reflections on the workspace using the given UB matrix
+ * @param peaksWS Workspace containing peaks
+ * @param ub A UB matrix to define the the transform from Q_sample to hkl
+ * @param tolerance If an index is within this tolerance of an integer then
+ * accept it
+ * @param optimizeUB If true optimize the UB for these peaks
+ * @param satelliteArgs If set, attempt to index peaks as satellies if main
+ * indexing fails
+ * @return A CombinedIndexingStats detailing the output found
+ */
+CombinedIndexingStats
+indexPeaks(const std::vector<Peak *> &peaks, DblMatrix ub,
+ const double mainTolerance, const bool roundHKLs,
+ const bool optimizeUB,
+ const Prop::SatelliteIndexingArgs &satelliteArgs) {
+ const auto nPeaks = peaks.size();
+ std::vector<V3D> qSample(nPeaks);
+ std::generate(
+ std::begin(qSample), std::end(qSample),
+ [&peaks, i = 0u ]() mutable { return peaks[i++]->getQSampleFrame(); });
+
+ if (optimizeUB) {
+ ub = optimizeUBMatrix(ub, qSample, mainTolerance);
+ }
- for (const auto &peak : peaks) {
- if (peak.getRunNumber() == run)
- q_vectors.push_back(peak.getQSampleFrame());
+ CombinedIndexingStats stats;
+ ub.Invert();
+ for (auto i = 0u; i < peaks.size(); ++i) {
+ const auto peak = peaks[i];
+ V3D millerIndices;
+ if (IndexingUtils::CalculateMillerIndices(ub, qSample[i], mainTolerance,
+ millerIndices)) {
+ stats.main.numIndexed++;
+ stats.main.error += millerIndices.hklError() / 3.0;
+ if (roundHKLs) {
+ IndexingUtils::RoundHKL(millerIndices);
}
+ peak->setHKL(millerIndices);
+ peak->setIntHKL(millerIndices);
+ peak->setIntMNP(V3D(0, 0, 0));
+ }
- Matrix<double> tempUB(UB);
+ // std::vector<V3D> millerIndices;
+ // millerIndices.reserve(qSample.size());
+ // CombinedIndexingStats stats;
+ // double averageError{0.0};
+ // stats.main.numIndexed = IndexingUtils::CalculateMillerIndices(
+ // ub, qSample, mainTolerance, millerIndices, averageError);
+ // stats.main.error = averageError * stats.main.numIndexed;
+
+ // if (static_cast<size_t>(stats.main.numIndexed) != nPeaks &&
+ // satelliteArgs.maxOrder > 0) {
+ // stats.satellites =
+ // indexSatellites(peaks, millerIndices, roundHKLs, satelliteArgs);
+ // } else {
+ // if (roundHKLs)
+ // IndexingUtils::RoundHKLs(millerIndices);
+
+ // std::for_each(std::begin(peaks), std::end(peaks),
+ // [&millerIndices, i = 0u ](Peak * peak) mutable {
+ // peak->setHKL(millerIndices[i]);
+ // peak->setIntHKL(millerIndices[i]);
+ // peak->setIntMNP(V3D(0, 0, 0));
+ // ++i;
+ // });
+ // }
+ }
+ return stats;
+}
- int num_indexed = 0;
- int original_indexed = 0;
- double original_error = 0;
- original_indexed = IndexingUtils::CalculateMillerIndices(
- tempUB, q_vectors, tolerance, miller_indices, original_error);
+} // namespace
- IndexingUtils::RoundHKLs(miller_indices); // HKLs must be rounded for
- // Optimize_UB to work
- num_indexed = original_indexed;
- average_error = original_error;
+/** Initialize the algorithm's properties.
+ */
+void IndexPeaks::init() {
+ using Mantid::API::WorkspaceProperty;
+ using Mantid::Kernel::BoundedValidator;
+ using Mantid::Kernel::Direction;
- bool done = false;
- if (num_indexed < 3) // can't optimize without at least 3
- { // peaks
- done = true;
- }
+ // -- inputs --
+ this->declareProperty(
+ std::make_unique<WorkspaceProperty<PeaksWorkspace_sptr::element_type>>(
+ Prop::PEAKSWORKSPACE, "", Direction::InOut),
+ "Input Peaks Workspace");
- int iteration = 0;
- while (iteration < 4 && !done) // try repeatedly optimizing 4 times
- { // which is usually sufficient
- try {
- IndexingUtils::Optimize_UB(tempUB, miller_indices, q_vectors);
- } catch (...) // If there is any problem, such as too few
- { // independent peaks, just use the original UB
- tempUB = UB;
- done = true;
- }
-
- num_indexed = IndexingUtils::CalculateMillerIndices(
- tempUB, q_vectors, tolerance, miller_indices, average_error);
-
- IndexingUtils::RoundHKLs(miller_indices); // HKLs must be rounded for
- // Optimize_UB to work
-
- if (num_indexed < original_indexed) // just use the original UB
- {
- num_indexed = original_indexed;
- average_error = original_error;
- done = true;
- }
-
- iteration++;
- }
+ auto mustBePositive = boost::make_shared<BoundedValidator<double>>();
+ mustBePositive->setLower(0.0);
+ this->declareProperty(Prop::TOLERANCE, 0.15, mustBePositive,
+ "Main peak indexing tolerance", Direction::Input);
+ this->declareProperty(Prop::SATE_TOLERANCE, 0.15, mustBePositive,
+ "Satellite peak indexing tolerance", Direction::Input);
+ this->declareProperty(Prop::ROUNDHKLS, true,
+ "Round H, K and L values to integers");
+ this->declareProperty(Prop::COMMONUB, false,
+ "Index all orientations with a common UB");
+ // -- outputs --
+ this->declareProperty(Prop::NUM_INDEXED, 0,
+ "Gets set with the number of indexed peaks.",
+ Direction::Output);
+ this->declareProperty(Prop::AVERAGE_ERR, 0.0,
+ "Gets set with the average HKL indexing error.",
+ Direction::Output);
+ this->declareProperty(Prop::MAIN_NUM_INDEXED, 0,
+ "Gets set with the number of indexed main peaks.",
+ Direction::Output);
+ this->declareProperty(Prop::SATE_NUM_INDEXED, 0,
+ "Gets set with the number of indexed main peaks.",
+ Direction::Output);
+ this->declareProperty(
+ Prop::MAIN_ERR, 0.0,
+ "Gets set with the average HKL indexing error of Main Peaks.",
+ Direction::Output);
+ this->declareProperty(
+ Prop::SATE_ERR, 0.0,
+ "Gets set with the average HKL indexing error of Satellite Peaks.",
+ Direction::Output);
+}
- // If data not modulated, recalculate fractional HKL
- if (o_lattice.getMaxOrder() == 0) {
- // If user wants fractional hkls, recalculate them
- if (!round_hkls) {
- num_indexed = IndexingUtils::CalculateMillerIndices(
- tempUB, q_vectors, tolerance, miller_indices, average_error);
- }
- total_indexed += num_indexed;
- total_error += average_error * num_indexed;
-
- // tell the user how many were indexed in each run
- if (run_numbers.size() > 1) {
- g_log.notice() << "Run " << run << ": indexed " << num_indexed
- << " Peaks out of " << q_vectors.size()
- << " with tolerance of " << tolerance << '\n';
- g_log.notice() << "Average error in h,k,l for indexed peaks = "
- << average_error << '\n';
- }
-
- size_t miller_index_counter = 0;
- for (auto &peak : peaks) {
- if (peak.getRunNumber() == run) {
- peak.setHKL(miller_indices[miller_index_counter]);
- peak.setIntHKL(miller_indices[miller_index_counter]);
- peak.setIntMNP(V3D(0, 0, 0));
- miller_index_counter++;
- }
- }
- } else {
- g_log.notice() << "Maximum Order: " << o_lattice.getMaxOrder() << '\n';
- int ModDim = 0;
- int main_indexed = 0;
- int sate_indexed = 0;
- double main_error = 0;
- double sate_error = 0;
- int maxOrder = o_lattice.getMaxOrder();
- bool crossTerm = o_lattice.getCrossTerm();
- V3D offsets1 = o_lattice.getModVec(0);
- V3D offsets2 = o_lattice.getModVec(1);
- V3D offsets3 = o_lattice.getModVec(2);
-
- if (offsets1 == V3D(0, 0, 0))
- throw std::runtime_error("Invalid Modulation Vector");
- else if (offsets2 == V3D(0, 0, 0))
- ModDim = 1;
- else if (offsets3 == V3D(0, 0, 0))
- ModDim = 2;
- else
- ModDim = 3;
-
- IndexingUtils::CalculateMillerIndices(tempUB, q_vectors, 1.0,
- miller_indices, average_error);
-
- // Index satellite peaks
- size_t miller_index_counter = 0;
- for (auto &peak : peaks) {
- if (peak.getRunNumber() == run) {
- peak.setHKL(miller_indices[miller_index_counter]);
- miller_index_counter++;
- auto hkl = peak.getHKL();
- bool peak_main_indexed{false}, peak_sat_indexed{false};
-
- if (IndexingUtils::ValidIndex(hkl, tolerance)) {
- if (round_hkls) {
- IndexingUtils::RoundHKL(hkl);
- peak.setHKL(hkl);
- }
- peak.setIntHKL(hkl);
- peak.setIntMNP(V3D(0, 0, 0));
- peak_main_indexed = true;
- main_error += hkl.hklError();
- } else if (!crossTerm) {
- if (ModDim > 0) {
- for (int order = -maxOrder; order <= maxOrder; order++) {
- if (order == 0)
- continue; // exclude order 0
- V3D hkl1(hkl);
- hkl1 -= offsets1 * order;
- if (IndexingUtils::ValidIndex(hkl1, satetolerance)) {
- peak.setIntHKL(hkl1);
- peak.setIntMNP(V3D(order, 0, 0));
- peak_sat_indexed = true;
- sate_error += hkl1.hklError();
- }
- }
- }
- if (ModDim > 1) {
- for (int order = -maxOrder; order <= maxOrder; order++) {
- if (order == 0)
- continue; // exclude order 0
- V3D hkl1(hkl);
- hkl1 -= offsets2 * order;
- if (IndexingUtils::ValidIndex(hkl1, satetolerance)) {
- peak.setIntHKL(hkl1);
- peak.setIntMNP(V3D(0, order, 0));
- peak_sat_indexed = true;
- sate_error += hkl1.hklError();
- }
- }
- }
- if (ModDim > 2) {
- for (int order = -maxOrder; order <= maxOrder; order++) {
- if (order == 0)
- continue; // exclude order 0
- V3D hkl1(hkl);
- hkl1 -= offsets3 * order;
- if (IndexingUtils::ValidIndex(hkl1, satetolerance)) {
- peak.setIntHKL(hkl1);
- peak.setIntMNP(V3D(0, 0, order));
- peak_sat_indexed = true;
- sate_error += hkl1.hklError();
- }
- }
- }
- } else {
- if (ModDim == 1) {
- for (int order = -maxOrder; order <= maxOrder; order++) {
- if (order == 0)
- continue; // exclude order 0
- V3D hkl1(hkl);
- hkl1 -= offsets1 * order;
- if (IndexingUtils::ValidIndex(hkl1, satetolerance)) {
- peak.setIntHKL(hkl1);
- peak.setIntMNP(V3D(order, 0, 0));
- peak_sat_indexed = true;
- sate_error += hkl1.hklError();
- }
- }
- }
- if (ModDim == 2) {
- for (int m = -maxOrder; m <= maxOrder; m++)
- for (int n = -maxOrder; n <= maxOrder; n++) {
- if (m == 0 && n == 0)
- continue; // exclude 0,0
- V3D hkl1(hkl);
- hkl1 -= offsets1 * m + offsets2 * n;
- if (IndexingUtils::ValidIndex(hkl1, satetolerance)) {
- peak.setIntHKL(hkl1);
- peak.setIntMNP(V3D(m, n, 0));
- peak_sat_indexed = true;
- sate_error += hkl1.hklError();
- }
- }
- }
- if (ModDim == 3) {
- for (int m = -maxOrder; m <= maxOrder; m++)
- for (int n = -maxOrder; n <= maxOrder; n++)
- for (int p = -maxOrder; p <= maxOrder; p++) {
- if (m == 0 && n == 0 && p == 0)
- continue; // exclude 0,0,0
- V3D hkl1(hkl);
- hkl1 -= offsets1 * m + offsets2 * n + offsets3 * p;
- if (IndexingUtils::ValidIndex(hkl1, satetolerance)) {
- peak.setIntHKL(hkl1);
- peak.setIntMNP(V3D(m, n, p));
- peak_sat_indexed = true;
- sate_error += hkl1.hklError();
- }
- }
- }
- }
- if (peak_main_indexed) {
- main_indexed++;
- peak.setIntHKL(V3D(0, 0, 0));
- peak.setIntMNP(V3D(0, 0, 0));
- } else if (peak_sat_indexed)
- sate_indexed++;
- else {
- peak.setHKL(V3D(0, 0, 0));
- peak.setIntHKL(V3D(0, 0, 0));
- peak.setIntMNP(V3D(0, 0, 0));
- }
- }
- }
-
- num_indexed = main_indexed + sate_indexed;
- total_main += main_indexed;
- total_sate += sate_indexed;
- total_main_error += main_error / 3;
- total_sate_error += sate_error / 3;
- total_indexed += main_indexed + sate_indexed;
- total_error += main_error / 3 + sate_error / 3;
-
- if (run_numbers.size() > 1) {
- g_log.notice() << "Run " << run << ": indexed " << num_indexed
- << " Peaks out of " << q_vectors.size() << '\n';
- g_log.notice() << "of which, " << main_indexed
- << " Main Bragg Peaks are indexed with tolerance of "
- << tolerance << ", " << sate_indexed
- << " Satellite Peaks are indexed with tolerance of "
- << satetolerance << '\n';
- }
- }
- }
+/**
+ * Validate all inputs once set
+ * @return A map of property name to help message
+ */
+std::map<std::string, std::string> IndexPeaks::validateInputs() {
+ std::map<std::string, std::string> helpMsgs;
+
+ PeaksWorkspace_sptr ws = this->getProperty(Prop::PEAKSWORKSPACE);
+ try {
+ ws->sample().getOrientedLattice();
+ } catch (std::runtime_error &exc) {
+ helpMsgs[Prop::PEAKSWORKSPACE] = exc.what();
+ }
- if (total_indexed > 0)
- average_error = total_error / total_indexed;
- else
- average_error = 0;
+ return helpMsgs;
+}
- if (total_main > 0)
- average_main_error = total_main_error / total_main;
- else
- average_main_error = 0;
+/**
+ * Execute the algorithm.
+ */
+void IndexPeaks::exec() {
+ const auto args = Prop::IndexPeaksArgs::parse(*this);
- if (total_sate > 0)
- average_sate_error = total_sate_error / total_sate;
- else
- average_sate_error = 0;
+ // quick exit
+ if (args.workspace->getNumberPeaks() == 0) {
+ g_log.warning("Empty peaks workspace. Nothing to index");
+ return;
}
- if (o_lattice.getMaxOrder() == 0) {
- // tell the user how many were indexed overall and the overall average error
- g_log.notice() << "ALL Runs: indexed " << total_indexed << " Peaks out of "
- << n_peaks << " with tolerance of " << tolerance << '\n';
- g_log.notice() << "Average error in h,k,l for indexed peaks = "
- << average_error << '\n';
-
- // Save output properties
- this->setProperty("NumIndexed", total_indexed);
- this->setProperty("MainNumIndexed", total_indexed);
- this->setProperty("SateNumIndexed", 0);
- this->setProperty("AverageError", average_error);
- this->setProperty("MainError", average_error);
- this->setProperty("SatelliteError", 0.);
- // Show the lattice parameters
- g_log.notice() << o_lattice << "\n";
+ CombinedIndexingStats indexingInfo;
+ const auto &lattice = args.workspace->sample().getOrientedLattice();
+ const auto &sampleUB = lattice.getUB();
+ if (args.commonUB) {
+ // Use sample UB an all peaks regardless of run
+ std::vector<Peak *> allPeaksRef(args.workspace->getNumberPeaks());
+ std::transform(std::begin(args.workspace->getPeaks()),
+ std::end(args.workspace->getPeaks()),
+ std::begin(allPeaksRef), [](Peak &peak) { return &peak; });
+ const bool optimizeUB{false};
+ indexingInfo = indexPeaks(allPeaksRef, sampleUB, args.mainTolerance,
+ args.roundHKLs, optimizeUB, args.satellites);
} else {
- g_log.notice() << "ALL Runs: indexed " << total_indexed << " Peaks out of "
- << n_peaks << " with tolerance of " << tolerance << '\n';
- g_log.notice() << "Out of " << total_indexed << " Indexed Peaks "
- << total_main << " are Main Bragg Peaks, and " << total_sate
+ // Use a UB optimized for each run
+ auto &allPeaks = args.workspace->getPeaks();
+ std::unordered_map<int, std::vector<Peak *>> peaksPerRun;
+ std::for_each(std::begin(allPeaks), std::end(allPeaks),
+ [&peaksPerRun](Peak &peak) {
+ peaksPerRun[peak.getRunNumber()].emplace_back(&peak);
+ });
+ const bool optimizeUB{true};
+ for (const auto &runPeaks : peaksPerRun) {
+ indexingInfo += indexPeaks(runPeaks.second, sampleUB, args.mainTolerance,
+ args.roundHKLs, optimizeUB, args.satellites);
+ }
+ }
+
+ // else {
+ // assert(false);
+ // auto ws = args.workspace;
+ // const auto &o_lattice = ws->sample().getOrientedLattice();
+ // const auto &UB = o_lattice.getUB();
+ // const bool round_hkls = this->getProperty(Prop::ROUNDHKLS);
+ // std::vector<Peak> &peaks = ws->getPeaks();
+ // const size_t n_peaks = static_cast<size_t>(ws->getNumberPeaks());
+ // int total_indexed{0}, total_main{0}, total_sate{0};
+ // double average_error{0.}, average_main_error{0.},
+ // average_sate_error{0.}; const double tolerance =
+ // this->getProperty("Tolerance");
+
+ // double total_error{0.}, total_main_error{0.}, total_sate_error{0.};
+ // double satetolerance = this->getProperty("ToleranceForSatellite");
+
+ // // get list of run numbers in this peaks workspace
+ // std::vector<int> run_numbers;
+ // for (const auto &peak : peaks) {
+ // int run = peak.getRunNumber();
+ // bool found = false;
+ // size_t k = 0;
+ // while (k < run_numbers.size() && !found) {
+ // if (run == run_numbers[k])
+ // found = true;
+ // else
+ // k++;
+ // }
+ // if (!found)
+ // run_numbers.push_back(run);
+ // }
+
+ // // index the peaks for each run separately, using a UB matrix
+ // optimized
+ // // for that run
+
+ // for (const int run : run_numbers) {
+ // std::vector<V3D> miller_indices;
+ // std::vector<V3D> q_vectors;
+
+ // for (const auto &peak : peaks) {
+ // if (peak.getRunNumber() == run)
+ // q_vectors.push_back(peak.getQSampleFrame());
+ // }
+
+ // DblMatrix tempUB(UB);
+
+ // int num_indexed = 0;
+ // double original_error = 0;
+ // int original_indexed = IndexingUtils::CalculateMillerIndices(
+ // tempUB, q_vectors, tolerance, miller_indices, original_error);
+
+ // IndexingUtils::RoundHKLs(miller_indices); // HKLs must be rounded
+ // for
+ // // Optimize_UB to work
+ // num_indexed = original_indexed;
+ // average_error = original_error;
+
+ // bool done = false;
+ // if (num_indexed < 3) // can't optimize without at least 3
+ // { // peaks
+ // done = true;
+ // }
+
+ // int iteration = 0;
+ // while (iteration < 4 && !done) // try repeatedly optimizing 4 times
+ // { // which is usually sufficient
+ // try {
+ // IndexingUtils::Optimize_UB(tempUB, miller_indices, q_vectors);
+ // } catch (...) // If there is any problem, such as too few
+ // { // independent peaks, just use the original UB
+ // tempUB = UB;
+ // done = true;
+ // }
+
+ // num_indexed = IndexingUtils::CalculateMillerIndices(
+ // tempUB, q_vectors, tolerance, miller_indices, average_error);
+
+ // IndexingUtils::RoundHKLs(miller_indices); // HKLs must be rounded
+ // for
+ // // Optimize_UB to work
+
+ // if (num_indexed < original_indexed) // just use the original UB
+ // {
+ // num_indexed = original_indexed;
+ // average_error = original_error;
+ // done = true;
+ // }
+
+ // iteration++;
+ // }
+
+ // // If data not modulated, recalculate fractional HKL
+ // if (o_lattice.getMaxOrder() == 0) {
+ // // If user wants fractional hkls, recalculate them
+ // if (!round_hkls) {
+ // num_indexed = IndexingUtils::CalculateMillerIndices(
+ // tempUB, q_vectors, tolerance, miller_indices,
+ // average_error);
+ // }
+ // total_indexed += num_indexed;
+ // total_main += num_indexed;
+ // total_error += average_error * num_indexed;
+ // total_main_error += average_error * num_indexed;
+
+ // // tell the user how many were indexed in each run
+ // if (run_numbers.size() > 1) {
+ // g_log.notice() << "Run " << run << ": indexed " << num_indexed
+ // << " Peaks out of " << q_vectors.size()
+ // << " with tolerance of " << tolerance << '\n';
+ // g_log.notice() << "Average error in h,k,l for indexed peaks = "
+ // << average_error << '\n';
+ // }
+
+ // size_t miller_index_counter = 0;
+ // for (auto &peak : peaks) {
+ // if (peak.getRunNumber() == run) {
+ // peak.setHKL(miller_indices[miller_index_counter]);
+ // peak.setIntHKL(miller_indices[miller_index_counter]);
+ // peak.setIntMNP(V3D(0, 0, 0));
+ // miller_index_counter++;
+ // }
+ // }
+
+ // } else {
+ // g_log.notice() << "Maximum Order: " << o_lattice.getMaxOrder() <<
+ // '\n'; int ModDim = 0; int main_indexed = 0; int sate_indexed = 0;
+ // double main_error = 0;
+ // double sate_error = 0;
+ // const int maxOrder = o_lattice.getMaxOrder();
+ // const bool crossTerm = o_lattice.getCrossTerm();
+ // const V3D offsets1 = o_lattice.getModVec(0);
+ // const V3D offsets2 = o_lattice.getModVec(1);
+ // const V3D offsets3 = o_lattice.getModVec(2);
+
+ // if (offsets1 == V3D(0, 0,
+ // 0))test_zero_satellite_tol_only_indexes_main_refl_with_modvectors_from_ub
+ // throw std::runtime_error("Invalid Modulation Vector");
+ // else if (offsets2 == V3D(0, 0, 0))
+ // ModDim = 1;
+ // else if (offsets3 == V3D(0, 0, 0))
+ // ModDim = 2;
+ // else
+ // ModDim = 3;
+
+ // IndexingUtils::CalculateMillerIndices(tempUB, q_vectors, 1.0,
+ // miller_indices,
+ // average_error);
+
+ // // Index satellite peaks
+ // size_t miller_index_counter = 0;
+ // for (auto &peak : peaks) {
+ // if (peak.getRunNumber() == run) {
+ // peak.setHKL(miller_indices[miller_index_counter]);
+ // miller_index_counter++;
+ // auto hkl = peak.getHKL();
+ // bool suc_indexed = false;
+
+ // if (IndexingUtils::ValidIndex(hkl, tolerance)) {
+ // if (round_hkls) {
+ // IndexingUtils::RoundHKL(hkl);
+ // peak.setHKL(hkl);
+ // }
+ // peak.setIntHKL(hkl);
+ // peak.setIntMNP(V3D(0, 0, 0));
+ // suc_indexed = true;
+ // main_indexed++;
+ // main_error += hkl.hklError();
+ // } else if (!crossTerm) {
+ // if (ModDim > 0) {
+ // for (int order = -maxOrder; order <= maxOrder; order++) {
+ // if (order == 0)
+ // continue; // exclude order 0
+ // V3D hkl1(hkl);
+ // hkl1 -= offsets1 * order;
+ // if (IndexingUtils::ValidIndex(hkl1, satetolerance)) {
+ // peak.setIntHKL(hkl1);
+ // peak.setIntMNP(V3D(order, 0, 0));
+ // suc_indexed = true;
+ // sate_indexed++;
+ // sate_error += hkl1.hklError();
+ // }
+ // }
+ // }
+ // if (ModDim > 1) {
+ // for (int order = -maxOrder; order <= maxOrder; order++) {
+ // if (order == 0)
+ // continue; // exclude order 0
+ // V3D hkl1(hkl);
+ // hkl1 -= offsets2 * order;
+ // if (IndexingUtils::ValidIndex(hkl1, satetolerance)) {
+ // peak.setIntHKL(hkl1);
+ // peak.setIntMNP(V3D(0, order, 0));
+ // suc_indexed = true;
+ // sate_indexed++;
+ // sate_error += hkl1.hklError();
+ // }
+ // }
+ // }
+ // if (ModDim > 2) {
+ // for (int order = -maxOrder; order <= maxOrder; order++) {
+ // if (order == 0)
+ // continue; // exclude order 0
+ // V3D hkl1(hkl);
+ // hkl1 -= offsets3 * order;
+ // if (IndexingUtils::ValidIndex(hkl1, satetolerance)) {
+ // peak.setIntHKL(hkl1);
+ // peak.setIntMNP(V3D(0, 0, order));
+ // suc_indexed = true;
+ // sate_indexed++;
+ // sate_error += hkl1.hklError();
+ // }
+ // }
+ // }
+ // } else {
+ // if (ModDim == 1) {
+ // for (int order = -maxOrder; order <= maxOrder; order++) {
+ // if (order == 0)
+ // continue; // exclude order 0
+ // V3D hkl1(hkl);
+ // hkl1 -= offsets1 * order;
+ // if (IndexingUtils::ValidIndex(hkl1, satetolerance)) {
+ // peak.setIntHKL(hkl1);
+ // peak.setIntMNP(V3D(order, 0, 0));
+ // suc_indexed = true;
+ // sate_indexed++;
+ // sate_error += hkl1.hklError();
+ // }
+ // }
+ // }
+ // if (ModDim == 2) {
+ // for (int m = -maxOrder; m <= maxOrder; m++)
+ // for (int n = -maxOrder; n <= maxOrder; n++) {
+ // if (m == 0 && n == 0)
+ // continue; // exclude 0,0
+ // V3D hkl1(hkl);
+ // hkl1 -= offsets1 * m + offsets2 * n;
+ // if (IndexingUtils::ValidIndex(hkl1, satetolerance)) {
+ // peak.setIntHKL(hkl1);
+ // peak.setIntMNP(V3D(m, n, 0));
+ // suc_indexed = true;
+ // sate_indexed++;
+ // sate_error += hkl1.hklError();
+ // }
+ // }
+ // }
+ // if (ModDim == 3) {
+ // for (int m = -maxOrder; m <= maxOrder; m++)
+ // for (int n = -maxOrder; n <= maxOrder; n++)
+ // for (int p = -maxOrder; p <= maxOrder; p++) {
+ // if (m == 0 && n == 0 && p == 0)
+ // continue; // exclude 0,0,0
+ // V3D hkl1(hkl);
+ // hkl1 -= offsets1 * m + offsets2 * n + offsets3 * p;
+ // if (IndexingUtils::ValidIndex(hkl1, satetolerance))
+ // {
+ // peak.setIntHKL(hkl1);
+ // peak.setIntMNP(V3D(m, n, p));
+ // suc_indexed = true;
+ // sate_indexed++;
+ // sate_error += hkl1.hklError();
+ // }
+ // }
+ // }
+ // }
+ // if (!suc_indexed) {
+ // peak.setHKL(V3D(0, 0, 0));
+ // peak.setIntHKL(V3D(0, 0, 0));
+ // peak.setIntMNP(V3D(0, 0, 0));
+ // }
+ // }
+ // }
+
+ // num_indexed = main_indexed + sate_indexed;
+ // total_main += main_indexed;
+ // total_sate += sate_indexed;
+ // total_main_error += main_error / 3;
+ // total_sate_error += sate_error / 3;
+ // total_indexed += main_indexed + sate_indexed;
+ // total_error += main_error / 3 + sate_error / 3;
+
+ // if (run_numbers.size() > 1) {
+ // g_log.notice() << "Run " << run << ": indexed " << num_indexed
+ // << " Peaks out of " << q_vectors.size() << '\n';
+ // g_log.notice() << "of which, " << main_indexed
+ // << " Main Bragg Peaks are indexed with tolerance
+ // of "
+ // << tolerance << ", " << sate_indexed
+ // << " Satellite Peaks are indexed with tolerance
+ // of "
+ // << satetolerance << '\n';
+ // }
+ // }
+ // }
+
+ // if (total_indexed > 0)
+ // average_error = total_error / total_indexed;
+ // else
+ // average_error = 0;
+
+ // if (total_main > 0)
+ // average_main_error = total_main_error / total_main;
+ // else
+ // average_main_error = 0;
+
+ // if (total_sate > 0)
+ // average_sate_error = total_sate_error / total_sate;
+ // else
+ // average_sate_error = 0;
+
+ // indexingInfo.numMainIndexed = total_main;
+ // indexingInfo.numSatIndexed = total_sate;
+ // indexingInfo.averageError = average_error;
+ // indexingInfo.mainError = average_main_error;
+ // indexingInfo.satelliteError = average_sate_error;
+ //}
+
+ setProperty("NumIndexed", indexingInfo.totalNumIndexed());
+ setProperty("MainNumIndexed", indexingInfo.main.numIndexed);
+ setProperty("SateNumIndexed", indexingInfo.satellites.numIndexed);
+ setProperty("AverageError", indexingInfo.averageError());
+ setProperty("MainError", indexingInfo.mainError());
+ setProperty("SatelliteError", indexingInfo.satelliteError());
+
+ if (args.satellites.maxOrder > 0) {
+ g_log.notice() << "ALL Runs: indexed " << indexingInfo.totalNumIndexed()
+ << " Peaks out of " << args.workspace->getNumberPeaks()
+ << " with tolerance of " << args.mainTolerance << '\n';
+ g_log.notice() << "Out of " << indexingInfo.totalNumIndexed()
+ << " Indexed Peaks " << indexingInfo.main.numIndexed
+ << " are Main Bragg Peaks, and "
+ << indexingInfo.satellites.numIndexed
<< " are satellite peaks " << '\n';
g_log.notice() << "Average error in h,k,l for indexed peaks = "
- << average_error << '\n';
+ << indexingInfo.averageError() << '\n';
g_log.notice() << "Average error in h,k,l for indexed main peaks = "
- << average_main_error << '\n';
+ << indexingInfo.main.error << '\n';
g_log.notice() << "Average error in h,k,l for indexed satellite peaks = "
- << average_sate_error << '\n';
-
- // Save output properties
- this->setProperty("NumIndexed", total_indexed);
- this->setProperty("MainNumIndexed", total_main);
- this->setProperty("SateNumIndexed", total_sate);
- this->setProperty("MainError", average_main_error);
- this->setProperty("SatelliteError", average_sate_error);
+ << indexingInfo.satellites.error << '\n';
+
+ // Show the lattice parameters
+ g_log.notice() << args.workspace->sample().getOrientedLattice() << "\n";
+ } else {
+ // tell the user how many were indexed overall and the overall average
+ // error
+ g_log.notice() << "ALL Runs: indexed " << indexingInfo.totalNumIndexed()
+ << " Peaks out of " << args.workspace->getNumberPeaks()
+ << " with tolerance of " << args.mainTolerance << '\n';
+ g_log.notice() << "Average error in h,k,l for indexed peaks = "
+ << indexingInfo.main.error << '\n';
+
// Show the lattice parameters
- g_log.notice() << o_lattice << "\n";
+ g_log.notice() << args.workspace->sample().getOrientedLattice() << "\n";
}
}
diff --git a/Framework/Crystal/test/IndexPeaksTest.h b/Framework/Crystal/test/IndexPeaksTest.h
index 4d943831b99..06321fb2711 100644
--- a/Framework/Crystal/test/IndexPeaksTest.h
+++ b/Framework/Crystal/test/IndexPeaksTest.h
@@ -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() {
diff --git a/Framework/DataObjects/inc/MantidDataObjects/Peak.h b/Framework/DataObjects/inc/MantidDataObjects/Peak.h
index 1852c99554d..cf8cfdaee31 100644
--- a/Framework/DataObjects/inc/MantidDataObjects/Peak.h
+++ b/Framework/DataObjects/inc/MantidDataObjects/Peak.h
@@ -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;
diff --git a/Framework/DataObjects/src/Peak.cpp b/Framework/DataObjects/src/Peak.cpp
index b9e2c5f54d3..12c9ad8ee8b 100644
--- a/Framework/DataObjects/src/Peak.cpp
+++ b/Framework/DataObjects/src/Peak.cpp
@@ -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; }
diff --git a/Framework/DataObjects/test/PeakTest.h b/Framework/DataObjects/test/PeakTest.h
index 4a879884234..0270964ba16 100644
--- a/Framework/DataObjects/test/PeakTest.h
+++ b/Framework/DataObjects/test/PeakTest.h
@@ -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);
diff --git a/Framework/Geometry/inc/MantidGeometry/Crystal/IPeak.h b/Framework/Geometry/inc/MantidGeometry/Crystal/IPeak.h
index 20486094132..7a602b5c5f4 100644
--- a/Framework/Geometry/inc/MantidGeometry/Crystal/IPeak.h
+++ b/Framework/Geometry/inc/MantidGeometry/Crystal/IPeak.h
@@ -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;
diff --git a/Framework/Geometry/inc/MantidGeometry/Crystal/IndexingUtils.h b/Framework/Geometry/inc/MantidGeometry/Crystal/IndexingUtils.h
index 6b46b339987..8f09e1e0c5f 100644
--- a/Framework/Geometry/inc/MantidGeometry/Crystal/IndexingUtils.h
+++ b/Framework/Geometry/inc/MantidGeometry/Crystal/IndexingUtils.h
@@ -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,
diff --git a/Framework/Geometry/src/Crystal/IndexingUtils.cpp b/Framework/Geometry/src/Crystal/IndexingUtils.cpp
index 1147400604b..6d07ac18d96 100644
--- a/Framework/Geometry/src/Crystal/IndexingUtils.cpp
+++ b/Framework/Geometry/src/Crystal/IndexingUtils.cpp
@@ -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.
*/
diff --git a/Framework/Geometry/test/IndexingUtilsTest.h b/Framework/Geometry/test/IndexingUtilsTest.h
index 35eabbf202e..535802c05b4 100644
--- a/Framework/Geometry/test/IndexingUtilsTest.h
+++ b/Framework/Geometry/test/IndexingUtilsTest.h
@@ -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};
diff --git a/Framework/Geometry/test/MockObjects.h b/Framework/Geometry/test/MockObjects.h
index 9df67d6b0ec..4aaaec4b1c9 100644
--- a/Framework/Geometry/test/MockObjects.h
+++ b/Framework/Geometry/test/MockObjects.h
@@ -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());
--
GitLab