From 5aedbdada32fbf34d95df958fcaf5c655cb75216 Mon Sep 17 00:00:00 2001
From: Ross Whitfield <whitfieldre@ornl.gov>
Date: Wed, 24 Mar 2021 11:06:54 -0400
Subject: [PATCH] Add unit and system tests for LeanElasticPeak in FindPeaksMD

---
 Framework/MDAlgorithms/test/FindPeaksMDTest.h | 107 ++++++++++++++++++
 .../framework/ConvertHFIRSCDtoMDETest.py      |  12 ++
 .../tests/framework/TOPAZPeakFinding.py       |  37 ++++++
 3 files changed, 156 insertions(+)

diff --git a/Framework/MDAlgorithms/test/FindPeaksMDTest.h b/Framework/MDAlgorithms/test/FindPeaksMDTest.h
index 0ba29c344a1..ac83974259c 100644
--- a/Framework/MDAlgorithms/test/FindPeaksMDTest.h
+++ b/Framework/MDAlgorithms/test/FindPeaksMDTest.h
@@ -7,6 +7,7 @@
 #pragma once
 
 #include "MantidAPI/FrameworkManager.h"
+#include "MantidDataObjects/LeanElasticPeaksWorkspace.h"
 #include "MantidDataObjects/PeaksWorkspace.h"
 #include "MantidKernel/PropertyWithValue.h"
 #include "MantidMDAlgorithms/FindPeaksMD.h"
@@ -247,6 +248,112 @@ public:
 
     AnalysisDataService::Instance().remove("MDEWS");
   }
+
+  void do_test_LeanElastic(bool expInfo, bool histo) {
+    FrameworkManager::Instance().exec(
+        "CreateMDWorkspace", 18, "Dimensions", "3", "EventType", "MDEvent",
+        "Extents", "-10,10,-10,10,-10,10", "Names",
+        "Q_sample_x,Q_sample_y,Q_sample_z", "Units", "-,-,-", "SplitInto", "5",
+        "SplitThreshold", "20", "MaxRecursionDepth", "15", "OutputWorkspace",
+        "MDEWS");
+
+    if (expInfo) {
+      Instrument_sptr inst =
+          ComponentCreationHelper::createTestInstrumentRectangular2(1, 100,
+                                                                    0.05);
+      IMDEventWorkspace_sptr ws;
+      TS_ASSERT_THROWS_NOTHING(
+          ws = AnalysisDataService::Instance().retrieveWS<IMDEventWorkspace>(
+              "MDEWS"));
+      ExperimentInfo_sptr ei(new ExperimentInfo());
+      ei->setInstrument(inst);
+      // Give it a run number
+      ei->mutableRun().addProperty(
+          new PropertyWithValue<std::string>("run_number", "12345"), true);
+      ws->addExperimentInfo(ei);
+    }
+
+    addPeak(1000, 1, 2, 3, 0.1);
+    addPeak(3000, 4, 5, 6, 0.2);
+    addPeak(5000, -5, -5, 5, 0.2);
+
+    if (histo) {
+      FrameworkManager::Instance().exec(
+          "BinMD", 14, "AxisAligned", "1", "AlignedDim0",
+          "Q_sample_x,-10,10,100", "AlignedDim1", "Q_sample_y,-10,10,100",
+          "AlignedDim2", "Q_sample_z,-10,10,100", "IterateEvents", "1",
+          "InputWorkspace", "MDEWS", "OutputWorkspace", "MDEWS");
+    }
+
+    std::string outWSName("peaksFound");
+    FindPeaksMD alg;
+    TS_ASSERT_THROWS_NOTHING(alg.initialize())
+    TS_ASSERT(alg.isInitialized())
+    TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("InputWorkspace", "MDEWS"));
+    TS_ASSERT_THROWS_NOTHING(
+        alg.setPropertyValue("OutputWorkspace", outWSName));
+    TS_ASSERT_THROWS_NOTHING(
+        alg.setPropertyValue("DensityThresholdFactor", "2.0"));
+    TS_ASSERT_THROWS_NOTHING(
+        alg.setPropertyValue("PeakDistanceThreshold", "0.7"));
+
+    if (expInfo)
+      TS_ASSERT_THROWS_NOTHING(
+          alg.setPropertyValue("OutputPeakType", "LeanElasticPeak"));
+
+    TS_ASSERT_THROWS_NOTHING(alg.execute(););
+    TS_ASSERT(alg.isExecuted());
+
+    // Retrieve the workspace from data service.
+    LeanElasticPeaksWorkspace_sptr ws;
+    TS_ASSERT_THROWS_NOTHING(
+        ws = AnalysisDataService::Instance()
+                 .retrieveWS<LeanElasticPeaksWorkspace>(outWSName));
+    TS_ASSERT(ws);
+    if (!ws)
+      return;
+
+    // Should find all 3 peaks.
+    TS_ASSERT_EQUALS(ws->getNumberPeaks(), 3);
+
+    TS_ASSERT_DELTA(ws->getPeak(0).getQSampleFrame()[0], -5.0, 0.11);
+    TS_ASSERT_DELTA(ws->getPeak(0).getQSampleFrame()[1], -5.0, 0.11);
+    TS_ASSERT_DELTA(ws->getPeak(0).getQSampleFrame()[2], 5.0, 0.11);
+    if (expInfo) {
+      TS_ASSERT_EQUALS(ws->getPeak(0).getRunNumber(), 12345);
+    } else {
+      TS_ASSERT_EQUALS(ws->getPeak(0).getRunNumber(), -1);
+    }
+    // Bin count = density of the box / 1e6
+    double BinCount = ws->getPeak(0).getBinCount();
+    if (histo) {
+      TS_ASSERT_DELTA(BinCount, 0.08375, 0.001);
+    } else {
+      TS_ASSERT_DELTA(BinCount, 7., 001000.);
+    }
+
+    TS_ASSERT_DELTA(ws->getPeak(1).getQSampleFrame()[0], 4.0, 0.11);
+    TS_ASSERT_DELTA(ws->getPeak(1).getQSampleFrame()[1], 5.0, 0.11);
+    TS_ASSERT_DELTA(ws->getPeak(1).getQSampleFrame()[2], 6.0, 0.11);
+
+    TS_ASSERT_DELTA(ws->getPeak(2).getQSampleFrame()[0], 1.0, 0.11);
+    TS_ASSERT_DELTA(ws->getPeak(2).getQSampleFrame()[1], 2.0, 0.11);
+    TS_ASSERT_DELTA(ws->getPeak(2).getQSampleFrame()[2], 3.0, 0.11);
+
+    AnalysisDataService::Instance().remove("MDEWS");
+  }
+
+  void test_exec_LeanElastic() { do_test_LeanElastic(false, false); }
+
+  void test_exec_LeanElastic_histo() { do_test_LeanElastic(false, true); }
+
+  void test_exec_LeanElastic_with_expInfo() {
+    do_test_LeanElastic(true, false);
+  }
+
+  void test_exec_LeanElastic_histo_with_expInfo() {
+    do_test_LeanElastic(true, true);
+  }
 };
 
 //=====================================================================================
diff --git a/Testing/SystemTests/tests/framework/ConvertHFIRSCDtoMDETest.py b/Testing/SystemTests/tests/framework/ConvertHFIRSCDtoMDETest.py
index d1bc1ab2b9c..0c278aadf23 100644
--- a/Testing/SystemTests/tests/framework/ConvertHFIRSCDtoMDETest.py
+++ b/Testing/SystemTests/tests/framework/ConvertHFIRSCDtoMDETest.py
@@ -45,6 +45,18 @@ class ConvertHFIRSCDtoMDETest(systemtesting.MantidSystemTest):
                                        atol=0.005,
                                        err_msg=f"mismatch for peak {p}")
 
+        # now try using LeanElasticPeak
+        ConvertHFIRSCDtoMDETest_peaks3 = FindPeaksMD(InputWorkspace=ConvertHFIRSCDtoMDETest_Q, PeakDistanceThreshold=2.2,
+                                                     OutputPeakType='LeanElasticPeak')
+
+        self.assertEqual(ConvertHFIRSCDtoMDETest_peaks3.getNumberPeaks(), 14)
+
+        for p in range(14):
+            np.testing.assert_allclose(ConvertHFIRSCDtoMDETest_peaks3.getPeak(p).getQSampleFrame(),
+                                       ConvertHFIRSCDtoMDETest_peaks.getPeak(p).getQSampleFrame(),
+                                       atol=0.005,
+                                       err_msg=f"mismatch for peak {p}")
+
 
 class ConvertHFIRSCDtoMDE_HB3A_Test(systemtesting.MantidSystemTest):
     def requiredMemoryMB(self):
diff --git a/Testing/SystemTests/tests/framework/TOPAZPeakFinding.py b/Testing/SystemTests/tests/framework/TOPAZPeakFinding.py
index 55e6b9354bf..89e7cc50c33 100644
--- a/Testing/SystemTests/tests/framework/TOPAZPeakFinding.py
+++ b/Testing/SystemTests/tests/framework/TOPAZPeakFinding.py
@@ -13,6 +13,7 @@ them.
 import systemtesting
 import numpy
 from mantid.simpleapi import *
+from mantid.dataobjects import PeaksWorkspace, LeanElasticPeaksWorkspace
 
 
 class TOPAZPeakFinding(systemtesting.MantidSystemTest):
@@ -72,6 +73,42 @@ class TOPAZPeakFinding(systemtesting.MantidSystemTest):
         ConvertToDiffractionMDWorkspace(InputWorkspace='topaz_3132',OutputWorkspace='topaz_3132_QSample',
                                         OutputDimensions='Q (sample frame)',LorentzCorrection='1',SplitInto='2',SplitThreshold='150')
         FindPeaksMD(InputWorkspace='topaz_3132_QSample',PeakDistanceThreshold='0.12',MaxPeaks='200',OutputWorkspace='peaks_QSample')
+        self.assertTrue(isinstance(mtd['peaks_QSample'], PeaksWorkspace))
+        FindUBUsingFFT(PeaksWorkspace='peaks_QSample',MinD='2',MaxD='16')
+        CopySample(InputWorkspace='peaks_QSample',OutputWorkspace='topaz_3132',CopyName='0',CopyMaterial='0',
+                   CopyEnvironment='0',CopyShape='0')
+
+        # Index the peaks and check
+        results = IndexPeaks(PeaksWorkspace='peaks_QSample')
+        indexed = results[0]
+        if indexed < 199:
+            raise Exception("Expected at least 199 of 200 peaks to be indexed. Only indexed %d!" % indexed)
+
+        # Check the UB matrix
+        w = mtd["topaz_3132"]
+        s = w.sample()
+        ol = s.getOrientedLattice()
+        self.assertDelta( ol.a(), 4.714, 0.01, "Correct lattice a value not found.")
+        self.assertDelta( ol.b(), 6.06, 0.01, "Correct lattice b value not found.")
+        self.assertDelta( ol.c(), 10.42, 0.01, "Correct lattice c value not found.")
+        self.assertDelta( ol.alpha(), 90, 0.4, "Correct lattice angle alpha value not found.")
+        self.assertDelta( ol.beta(), 90, 0.4, "Correct lattice angle beta value not found.")
+        self.assertDelta( ol.gamma(), 90, 0.4, "Correct lattice angle gamma value not found.")
+
+        # Compare new and old UBs
+        newUB = numpy.array(mtd["topaz_3132"].sample().getOrientedLattice().getUB())
+        # UB Matrices are not necessarily the same, some of the H,K and/or L sign can be reversed
+        diff = abs(newUB) - abs(originalUB) < 0.001
+        for c in range(3):
+            # This compares each column, allowing old == new OR old == -new
+            if not numpy.all(diff[:,c]) :
+                raise Exception("More than 0.001 difference between UB matrices: Q (lab frame):\n"
+                                "%s\nQ (sample frame):\n%s" % (originalUB, newUB) )
+
+        # repeat but use LeanElasticPeaks
+        FindPeaksMD(InputWorkspace='topaz_3132_QSample',PeakDistanceThreshold='0.12',MaxPeaks='200',OutputWorkspace='peaks_QSample',
+                    OutputPeakType='LeanElasticPeak')
+        self.assertTrue(isinstance(mtd['peaks_QSample'], LeanElasticPeaksWorkspace))
         FindUBUsingFFT(PeaksWorkspace='peaks_QSample',MinD='2',MaxD='16')
         CopySample(InputWorkspace='peaks_QSample',OutputWorkspace='topaz_3132',CopyName='0',CopyMaterial='0',
                    CopyEnvironment='0',CopyShape='0')
-- 
GitLab