diff --git a/Framework/Crystal/src/IndexPeaks.cpp b/Framework/Crystal/src/IndexPeaks.cpp index dd851aadb089a9a9a463d4768581572cd69bd0f7..6cf2e353ab93ca7960308defbd48df95eafac2b8 100644 --- a/Framework/Crystal/src/IndexPeaks.cpp +++ b/Framework/Crystal/src/IndexPeaks.cpp @@ -47,17 +47,13 @@ void IndexPeaks::init() { this->declareProperty("CommonUBForAll", false, "Index all orientations with a common UB"); - this->declareProperty(std::make_unique<PropertyWithValue<int>>( - "NumIndexed", 0, Direction::Output), - "Gets set with the number of indexed peaks."); - this->declareProperty(std::make_unique<PropertyWithValue<double>>( "AverageError", 0.0, Direction::Output), "Gets set with the average HKL indexing error."); this->declareProperty(std::make_unique<PropertyWithValue<int>>( - "TotalNumIndexed", 0, Direction::Output), - "Gets set with the number of Total indexed peaks."); + "NumIndexed", 0, Direction::Output), + "Gets set with the number of indexed peaks."); this->declareProperty(std::make_unique<PropertyWithValue<int>>( "MainNumIndexed", 0, Direction::Output), @@ -82,10 +78,6 @@ void IndexPeaks::init() { */ void IndexPeaks::exec() { PeaksWorkspace_sptr ws = this->getProperty("PeaksWorkspace"); - if (!ws) { - throw std::runtime_error("Could not read the peaks workspace"); - } - OrientedLattice o_lattice = ws->mutableSample().getOrientedLattice(); const Matrix<double> &UB = o_lattice.getUB(); @@ -118,6 +110,8 @@ void IndexPeaks::exec() { 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]); @@ -204,11 +198,10 @@ void IndexPeaks::exec() { iteration++; } - if (o_lattice.getMaxOrder() == - 0) // If data not modulated, recalculate fractional HKL - { - if (!round_hkls) // If user wants fractional hkls, recalculate them - { + // 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); } @@ -264,11 +257,14 @@ void IndexPeaks::exec() { 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; @@ -422,7 +418,11 @@ void IndexPeaks::exec() { // 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"; } else { @@ -440,11 +440,11 @@ void IndexPeaks::exec() { << average_sate_error << '\n'; // Save output properties - setProperty("TotalNumIndexed", total_indexed); - setProperty("MainNumIndexed", total_main); - setProperty("SateNumIndexed", total_sate); - setProperty("MainError", average_main_error); - setProperty("SatelliteError", average_sate_error); + 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); // Show the lattice parameters g_log.notice() << o_lattice << "\n"; } diff --git a/Framework/Crystal/test/IndexPeaksTest.h b/Framework/Crystal/test/IndexPeaksTest.h index 9ff2d2b864bc32e3e18fc33c155556fd4dd50636..4d943831b994e684fdb6ca2facfbd607a60ce856 100644 --- a/Framework/Crystal/test/IndexPeaksTest.h +++ b/Framework/Crystal/test/IndexPeaksTest.h @@ -7,25 +7,109 @@ #ifndef MANTID_CRYSTAL_INDEX_PEAKS_TEST_H_ #define MANTID_CRYSTAL_INDEX_PEAKS_TEST_H_ -#include "MantidAPI/AnalysisDataService.h" #include "MantidAPI/Sample.h" -#include "MantidKernel/System.h" -#include "MantidKernel/Timer.h" -#include <cxxtest/TestSuite.h> - #include "MantidCrystal/IndexPeaks.h" -#include "MantidCrystal/LoadIsawUB.h" -#include "MantidDataHandling/LoadNexusProcessed.h" #include "MantidDataObjects/PeaksWorkspace.h" #include "MantidGeometry/Crystal/IndexingUtils.h" #include "MantidGeometry/Crystal/OrientedLattice.h" +#include "MantidTestHelpers/WorkspaceCreationHelper.h" + +#include <cxxtest/TestSuite.h> + +using Mantid::API::Workspace_sptr; +using Mantid::Crystal::IndexPeaks; +using Mantid::DataObjects::PeaksWorkspace_sptr; +using Mantid::Geometry::IndexingUtils; +using Mantid::Geometry::OrientedLattice; +using Mantid::Kernel::V3D; + +namespace { +// run number, det pos, q_sample +using MinimalPeak = std::tuple<int, V3D>; +template <int NPeaks> using MinimalPeaksList = std::array<MinimalPeak, NPeaks>; + +template <int NPeaks> +PeaksWorkspace_sptr +createPeaksWorkspace(const MinimalPeaksList<NPeaks> &testPeaksInfo, + const std::vector<double> &ub, const int maxOrder = -1, + const std::vector<V3D> &modVectors = std::vector<V3D>()) { + auto peaksWS = WorkspaceCreationHelper::createPeaksWorkspace(NPeaks); + OrientedLattice lattice; + lattice.setUB(ub); + if (maxOrder > 0) { + lattice.setMaxOrder(maxOrder); + if (modVectors.empty() || modVectors.size() > 3) + throw std::runtime_error("Only <= 3 modulation vectors can be provided."); + auto modVecOrDefault = [&modVectors](const size_t index) { + if (index < modVectors.size()) + return modVectors[index]; + else + return V3D(); + }; + lattice.setModVec1(modVecOrDefault(0)); + lattice.setModVec2(modVecOrDefault(1)); + lattice.setModVec3(modVecOrDefault(2)); + } + peaksWS->mutableSample().setOrientedLattice(&lattice); + + int peakCount{0}; + for (const auto &testPeakInfo : testPeaksInfo) { + auto &wsPeak = peaksWS->getPeak(peakCount); + wsPeak.setRunNumber(std::get<0>(testPeakInfo)); + wsPeak.setQSampleFrame(std::get<1>(testPeakInfo)); + ++peakCount; + } + + return peaksWS; +} + +PeaksWorkspace_sptr createTestPeaksWorkspaceMainReflOnly() { + constexpr int npeaks{5}; + const std::vector<double> ub = {-0.0122354, 0.00480056, -0.0860404, + 0.1165450, 0.00178145, 0.0045884, + 0.0273738, -0.08973560, 0.0252595}; + // peaks from TOPAZ_3007.peaks: 0, 1, 2, 10, 42 with sign for Q swapped + // as we don't use the crystallographic convention + constexpr std::array<MinimalPeak, npeaks> testPeaksInfo = { + MinimalPeak{3007, V3D(-3.52961, 3.13589, 1.0899)}, + MinimalPeak{3007, V3D(-2.42456, 2.29581, 1.71147)}, + MinimalPeak{3007, V3D(-3.04393, 3.05739, 2.03727)}, + MinimalPeak{3007, V3D(-4.02271, 2.4073, 1.62228)}, + MinimalPeak{3007, V3D(-4.04552, 1.59916, 3.71776)}}; + return createPeaksWorkspace<npeaks>(testPeaksInfo, ub); +} + +PeaksWorkspace_sptr +createTestPeaksWorkspaceWithSatellites(const int maxOrder, + const std::vector<V3D> &modVectors) { + constexpr int npeaks{5}; + const std::vector<double> ub = {0.269, -0.01, 0.033, 0.081, -0.191, + -0.039, 0.279, 0.347, -0.02}; + constexpr std::array<MinimalPeak, npeaks> testPeaksInfo = { + MinimalPeak{1, V3D(-3.691, -0.694, 3.762)}, // main + MinimalPeak{1, V3D(-1.234, -0.225, 1.25212)}, // satellite + MinimalPeak{1, V3D(-3.824, 0.728, 1.711)}, // main + MinimalPeak{1, V3D(0.872, -0.1998, 2.7476)}, // satellite + MinimalPeak{1, V3D(-1.54093, 0.129343, 1.445)}, // satellite + }; + return createPeaksWorkspace<npeaks>(testPeaksInfo, ub, maxOrder, modVectors); +} + +std::unique_ptr<IndexPeaks> +indexPeaks(PeaksWorkspace_sptr peaksWS, + const std::unordered_map<std::string, std::string> &arguments) { + auto alg = std::make_unique<IndexPeaks>(); + alg->setChild(true); + alg->initialize(); + alg->setProperty("PeaksWorkspace", peaksWS); + for (const auto &argument : arguments) { + alg->setPropertyValue(argument.first, argument.second); + } + alg->execute(); + return alg; +} -using namespace Mantid::Crystal; -using namespace Mantid::API; -using namespace Mantid::DataObjects; -using Mantid::DataHandling::LoadNexusProcessed; -using namespace Mantid::Kernel; -using namespace Mantid::Geometry; +} // namespace class IndexPeaksTest : public CxxTest::TestSuite { public: @@ -35,143 +119,205 @@ public: TS_ASSERT(alg.isInitialized()) } - void test_exec() { - // Name of the output workspace. - std::string WSName("peaks"); - LoadNexusProcessed loader; - TS_ASSERT_THROWS_NOTHING(loader.initialize()); - TS_ASSERT(loader.isInitialized()); - loader.setPropertyValue("Filename", "TOPAZ_3007.peaks.nxs"); - 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); - // clear all the peak indexes - std::vector<Peak> &peaks = ws->getPeaks(); - size_t n_peaks = ws->getNumberPeaks(); - for (size_t i = 0; i < n_peaks; i++) { - peaks[i].setHKL(V3D(0, 0, 0)); - } - // now set a proper UB in the - // oriented lattice - V3D row_0(-0.0122354, 0.00480056, -0.0860404); - V3D row_1(0.1165450, 0.00178145, 0.0045884); - V3D row_2(0.0273738, -0.08973560, 0.0252595); - - Matrix<double> UB(3, 3, false); - UB.setRow(0, row_0); - UB.setRow(1, row_1); - UB.setRow(2, row_2); - - OrientedLattice o_lattice; - o_lattice.setUB(UB); - ws->mutableSample().setOrientedLattice(&o_lattice); - - // index the peaks with the new UB - IndexPeaks alg; - TS_ASSERT_THROWS_NOTHING(alg.initialize()) - TS_ASSERT(alg.isInitialized()) - TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("PeaksWorkspace", WSName)); - TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("Tolerance", "0.1")); - TS_ASSERT_THROWS_NOTHING(alg.setProperty("RoundHKLs", false)); - TS_ASSERT_THROWS_NOTHING(alg.execute();); - TS_ASSERT(alg.isExecuted()); + void test_empty_peaks_workspace_indexes_nothing() { + auto peaksWS = WorkspaceCreationHelper::createPeaksWorkspace(0); + OrientedLattice lattice; + peaksWS->mutableSample().setOrientedLattice(&lattice); + + auto alg = indexPeaks(peaksWS, {{"Tolerance", "0.1"}}); - double tolerance = alg.getProperty("Tolerance"); + TS_ASSERT(alg->isExecuted()) + assertNumberPeaksIndexed(*alg, 0, 0, 0); + } + void test_no_roundHKLS_leaves_peaks_as_found() { + const auto ws = createTestPeaksWorkspaceMainReflOnly(); + const auto mainToleranceAsStr = "0.1"; + auto alg = + indexPeaks(ws, {{"Tolerance", mainToleranceAsStr}, {"RoundHKLs", "0"}}); // Check that the peaks were all indexed - for (size_t i = 0; i < n_peaks; i++) { - TS_ASSERT(IndexingUtils::ValidIndex(peaks[i].getHKL(), tolerance)); + const double mainTolerance{std::stod(mainToleranceAsStr)}; + auto &peaks = ws->getPeaks(); + for (const auto &peak : peaks) { + TS_ASSERT(IndexingUtils::ValidIndex(peak.getHKL(), mainTolerance)) } // Check the output properties - int numIndexed = alg.getProperty("NumIndexed"); - TS_ASSERT_EQUALS(numIndexed, 43); + assertNumberPeaksIndexed(*alg, 5, 5, 0); - double averageError = alg.getProperty("AverageError"); - TS_ASSERT_DELTA(averageError, 0.0097, 1e-3); + const double averageError = alg->getProperty("AverageError"); + TS_ASSERT_DELTA(averageError, 0.00639, 1e-4) // spot check a few peaks for // fractional Miller indices - V3D peak_0_hkl_d(-4.03065, -0.9885090, -6.01095); // first peak - V3D peak_1_hkl_d(-2.99276, 0.9955220, -4.00375); - V3D peak_2_hkl_d(-3.99311, 0.9856010, -5.00772); - V3D peak_10_hkl_d(-3.01107, -0.0155531, -7.01377); - V3D peak_42_hkl_d(-1.97065, 4.0283600, -6.97828); // last peak + 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 V3D error = peak_0_hkl_d - peaks[0].getHKL(); - TS_ASSERT_DELTA(error.norm(), 0.0, 2e-4); + TS_ASSERT_DELTA(error.norm(), 0.0, 2e-4) error = peak_1_hkl_d - peaks[1].getHKL(); - TS_ASSERT_DELTA(error.norm(), 0.0, 1e-4); + TS_ASSERT_DELTA(error.norm(), 0.0, 1e-4) error = peak_2_hkl_d - peaks[2].getHKL(); - TS_ASSERT_DELTA(error.norm(), 0.0, 1e-4); + TS_ASSERT_DELTA(error.norm(), 0.0, 1e-4) - error = peak_10_hkl_d - peaks[10].getHKL(); - TS_ASSERT_DELTA(error.norm(), 0.0, 1e-4); + error = peak_3_hkl_d - peaks[3].getHKL(); + TS_ASSERT_DELTA(error.norm(), 0.0, 1e-4) - error = peak_42_hkl_d - peaks[42].getHKL(); - TS_ASSERT_DELTA(error.norm(), 0.0, 2e-4); + error = peak_4_hkl_d - peaks[4].getHKL(); + TS_ASSERT_DELTA(error.norm(), 0.0, 2e-4) + } - // clear all the peak indexes then - // re-run the algorithm, rounding the HKLs - // this time, and again check a few peaks - for (size_t i = 0; i < n_peaks; i++) { - peaks[i].setHKL(V3D(0, 0, 0)); + void test_roundHKLS_gives_integer_miller_indices() { + const auto ws = createTestPeaksWorkspaceMainReflOnly(); + const auto mainToleranceAsStr = "0.1"; + auto alg = + indexPeaks(ws, {{"Tolerance", mainToleranceAsStr}, {"RoundHKLs", "1"}}); + // 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)) } - TS_ASSERT_THROWS_NOTHING(alg.initialize()) - TS_ASSERT(alg.isInitialized()) - TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("PeaksWorkspace", WSName)); - TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("Tolerance", "0.1")); - TS_ASSERT_THROWS_NOTHING(alg.setProperty("RoundHKLs", true)); - TS_ASSERT_THROWS_NOTHING(alg.execute();); - TS_ASSERT(alg.isExecuted()); - // Check that the peaks were all indexed - for (size_t i = 0; i < n_peaks; i++) { - TS_ASSERT(IndexingUtils::ValidIndex(peaks[i].getHKL(), tolerance)); + for (const auto &peak : peaks) { + TS_ASSERT(IndexingUtils::ValidIndex(peak.getHKL(), mainTolerance)) } // Check the output properties - numIndexed = alg.getProperty("NumIndexed"); - TS_ASSERT_EQUALS(numIndexed, 43); + assertNumberPeaksIndexed(*alg, 5, 5, 0); - averageError = alg.getProperty("AverageError"); - TS_ASSERT_DELTA(averageError, 0.0097, 1e-3); + const double averageError = alg->getProperty("AverageError"); + TS_ASSERT_DELTA(averageError, 0.00673, 1e-3) // spot check a few peaks for // integer 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_10_hkl(-3, 0, -7); - V3D peak_42_hkl(-2, 4, -7); // last peak + 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); // last peak - error = peak_0_hkl - peaks[0].getHKL(); - TS_ASSERT_DELTA(error.norm(), 0.0, 1e-10); + auto error = peak_0_hkl - peaks[0].getHKL(); + TS_ASSERT_DELTA(error.norm(), 0.0, 1e-10) error = peak_1_hkl - peaks[1].getHKL(); - TS_ASSERT_DELTA(error.norm(), 0.0, 1e-10); + TS_ASSERT_DELTA(error.norm(), 0.0, 1e-10) error = peak_2_hkl - peaks[2].getHKL(); - TS_ASSERT_DELTA(error.norm(), 0.0, 1e-10); + TS_ASSERT_DELTA(error.norm(), 0.0, 1e-10) + + error = peak_3_hkl - peaks[3].getHKL(); + TS_ASSERT_DELTA(error.norm(), 0.0, 1e-10) + + error = peak_4_hkl - peaks[4].getHKL(); + TS_ASSERT_DELTA(error.norm(), 0.0, 1e-10) + } - error = peak_10_hkl - peaks[10].getHKL(); - TS_ASSERT_DELTA(error.norm(), 0.0, 1e-10); + void + test_zero_satellite_tol_only_indexes_main_refl_with_modvectors_on_from_ub() { + const auto peaksWS = + createTestPeaksWorkspaceWithSatellites(1, {V3D(0.333, -0.667, 0.333)}); + const auto mainTolerance{"0.1"}, sateTolerance{"0."}; - error = peak_42_hkl - peaks[42].getHKL(); - TS_ASSERT_DELTA(error.norm(), 0.0, 1e-10); + const auto alg = + indexPeaks(peaksWS, {{"Tolerance", mainTolerance}, + {"ToleranceForSatellite", sateTolerance}}); + + assertNumberPeaksIndexed(*alg, 2, 2, 0); + assertIndexesAsExpected(*peaksWS, + {V3D(-1, 2, -9), V3D(0, 0, 0), V3D(-1, 1, -10), + V3D(0, 0, 0), V3D(0, 0, 0)}); + } + + void test_zero_maxorder_on_ub_indexes_main_refl_only() { + const auto peaksWS = + createTestPeaksWorkspaceWithSatellites(0, std::vector<V3D>()); + const auto mainTolerance{"0.1"}, sateTolerance{"1."}; + + const auto alg = + indexPeaks(peaksWS, {{"Tolerance", mainTolerance}, + {"ToleranceForSatellite", sateTolerance}}); + + assertNumberPeaksIndexed(*alg, 2, 2, 0); + assertIndexesAsExpected(*peaksWS, + {V3D(-1, 2, -9), V3D(0, 0, 0), V3D(-1, 1, -10), + V3D(0, 0, 0), V3D(0, 0, 0)}); + } + + void test_exec_with_common_ub_indexes_main_reflections_only() { + const auto peaksWS = + createTestPeaksWorkspaceWithSatellites(0, std::vector<V3D>()); + const auto sateTolerance{"1."}; + + const auto alg = + indexPeaks(peaksWS, {{"CommonUBForAll", "1"}, + {"ToleranceForSatellite", sateTolerance}}); + + assertNumberPeaksIndexed(*alg, 2, 2, 0); + 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() { + const auto peaksWS = WorkspaceCreationHelper::createPeaksWorkspace(1); + IndexPeaks alg; + alg.initialize(); + alg.setProperty("PeaksWorkspace", peaksWS); + + auto helpMsgs = alg.validateInputs(); + + const auto valueIter = helpMsgs.find("PeaksWorkspace"); + TS_ASSERT(valueIter != helpMsgs.cend()) + if (valueIter != helpMsgs.cend()) { + const auto msg = valueIter->second; + TS_ASSERT(msg.find("OrientedLattice") != std::string::npos) + } + } + +private: + void assertIndexesAsExpected(const PeaksWorkspace_sptr::element_type &peaksWS, + const std::vector<V3D> &expectedIndexes) { + assert(static_cast<size_t>(peaksWS.getNumberPeaks()) == + expectedIndexes.size()); + for (auto i = 0u; i < expectedIndexes.size(); ++i) { + TS_ASSERT_EQUALS(expectedIndexes[i], peaksWS.getPeak(i).getHKL()) + } + } + + void assertNumberPeaksIndexed(const IndexPeaks &alg, + const int numIndexedExpected, + const int mainNumIndexedExpected, + const int sateNumIndexedExpected) { + TS_ASSERT_EQUALS(numIndexedExpected, + static_cast<int>(alg.getProperty("NumIndexed"))) + TS_ASSERT_EQUALS(mainNumIndexedExpected, + static_cast<int>(alg.getProperty("MainNumIndexed"))) + TS_ASSERT_EQUALS(sateNumIndexedExpected, + static_cast<int>(alg.getProperty("SateNumIndexed"))) + } - // Remove workspace from the data service. - AnalysisDataService::Instance().remove(WSName); + void assertErrorsAsExpected(const IndexPeaks &alg, + const double expectedAverageError, + const double expectedMainError, + const double expectedSatError) { + TS_ASSERT_DELTA(expectedAverageError, + static_cast<double>(alg.getProperty("AverageError")), 1e-5); + TS_ASSERT_DELTA(expectedMainError, + static_cast<double>(alg.getProperty("MainError")), 1e-5); + TS_ASSERT_DELTA(expectedSatError, + static_cast<double>(alg.getProperty("SatelliteError")), + 1e-5); } };