Merge pull request #30953 from rosswhitfield/FindPeaksMD_LeanElasticPeaks

Modify FindPeaksMD to optionally produce LeanElasticPeaksWorkspace

Framework/MDAlgorithms/inc/MantidMDAlgorithms/FindPeaksMD.h

@@ -9,7 +9,9 @@
 #include "MantidAPI/Algorithm.h"
 #include "MantidAPI/ExperimentInfo.h"
 #include "MantidAPI/IMDEventWorkspace_fwd.h"
+#include "MantidAPI/IPeaksWorkspace.h"
 #include "MantidAPI/Progress.h"
+#include "MantidDataObjects/LeanElasticPeaksWorkspace.h"
 #include "MantidDataObjects/MDEventWorkspace.h"
 #include "MantidDataObjects/MDHistoWorkspace.h"
 #include "MantidDataObjects/PeaksWorkspace.h"
@@ -58,18 +60,29 @@ private:
   void exec() override;
 
   /// Read member variables from experiment info
-  void readExperimentInfo(const Mantid::API::ExperimentInfo_sptr &ei,
-                          const Mantid::API::IMDWorkspace_sptr &ws);
+  void readExperimentInfo(const Mantid::API::ExperimentInfo_sptr &ei);
+  void checkWorkspaceDims(const Mantid::API::IMDWorkspace_sptr &ws);
+  void determineOutputType(const std::string peakType,
+                           const uint16_t numExperimentInfo);
 
   /// Adds a peak based on Q, bin count & a set of detector IDs
   void addPeak(const Mantid::Kernel::V3D &Q, const double binCount,
                const Geometry::InstrumentRayTracer &tracer);
 
+  /// Adds a peak based on Q, bin count
+  void addLeanElasticPeak(const Mantid::Kernel::V3D &Q, const double binCount,
+                          const bool useGoniometer = false);
+
   /// Adds a peak based on Q, bin count
   std::shared_ptr<DataObjects::Peak>
   createPeak(const Mantid::Kernel::V3D &Q, const double binCount,
              const Geometry::InstrumentRayTracer &tracer);
 
+  /// Adds a peak based on Q, bin count
+  std::shared_ptr<DataObjects::LeanElasticPeak>
+  createLeanElasticPeak(const Mantid::Kernel::V3D &Q, const double binCount,
+                        const bool useGoniometer = false);
+
   /// Run find peaks on an MDEventWorkspace
   template <typename MDE, size_t nd>
   void findPeaks(typename DataObjects::MDEventWorkspace<MDE, nd>::sptr ws);
@@ -77,7 +90,7 @@ private:
   void findPeaksHisto(const Mantid::DataObjects::MDHistoWorkspace_sptr &ws);
 
   /// Output PeaksWorkspace
-  Mantid::DataObjects::PeaksWorkspace_sptr peakWS;
+  Mantid::API::IPeaksWorkspace_sptr peakWS;
 
   /// Estimated radius of peaks. Boxes closer than this are rejected
   coord_t peakRadiusSquared;
@@ -121,6 +134,8 @@ private:
   static const std::string volumeNormalization;
   /// NumberOfEventNormalization
   static const std::string numberOfEventsNormalization;
+
+  bool m_leanElasticPeak = false;
 };
 
 } // namespace MDAlgorithms

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"
@@ -29,7 +30,7 @@ static void createMDEW() {
       "CreateMDWorkspace", 18, "Dimensions", "3", "EventType", "MDEvent",
       "Extents", "-10,10,-10,10,-10,10", "Names", "Q_lab_x,Q_lab_y,Q_lab_z",
       "Units", "-,-,-", "SplitInto", "5", "SplitThreshold", "20",
-      "MaxRecursionDepth", "15", "OutputWorkspace", "MDEWS");
+      "MaxRecursionDepth", "15", "OutputWorkspace", "MDWS");
 
   // Give it an instrument
   Instrument_sptr inst =
@@ -37,7 +38,7 @@ static void createMDEW() {
   IMDEventWorkspace_sptr ws;
   TS_ASSERT_THROWS_NOTHING(
       ws = AnalysisDataService::Instance().retrieveWS<IMDEventWorkspace>(
-          "MDEWS"));
+          "MDWS"));
   ExperimentInfo_sptr ei(new ExperimentInfo());
   ei->setInstrument(inst);
   // Give it a run number
@@ -52,14 +53,14 @@ static void addPeak(size_t num, double x, double y, double z, double radius) {
   std::ostringstream mess;
   mess << num / 2 << ", " << x << ", " << y << ", " << z << ", " << radius;
   FrameworkManager::Instance().exec("FakeMDEventData", 4, "InputWorkspace",
-                                    "MDEWS", "PeakParams", mess.str().c_str());
+                                    "MDWS", "PeakParams", mess.str().c_str());
 
   // Add a center with more events (half radius, half the total), to create a
   // "peak"
   std::ostringstream mess2;
   mess2 << num / 2 << ", " << x << ", " << y << ", " << z << ", " << radius / 2;
   FrameworkManager::Instance().exec("FakeMDEventData", 4, "InputWorkspace",
-                                    "MDEWS", "PeakParams", mess2.str().c_str());
+                                    "MDWS", "PeakParams", mess2.str().c_str());
 }
 
 //=====================================================================================
@@ -91,14 +92,14 @@ public:
       FrameworkManager::Instance().exec(
           "BinMD", 14, "AxisAligned", "1", "AlignedDim0", "Q_lab_x,-10,10,100",
           "AlignedDim1", "Q_lab_y,-10,10,100", "AlignedDim2",
-          "Q_lab_z,-10,10,100", "IterateEvents", "1", "InputWorkspace", "MDEWS",
-          "OutputWorkspace", "MDEWS");
+          "Q_lab_z,-10,10,100", "IterateEvents", "1", "InputWorkspace", "MDWS",
+          "OutputWorkspace", "MDWS");
     }
 
     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("InputWorkspace", "MDWS"));
     TS_ASSERT_THROWS_NOTHING(
         alg.setPropertyValue("OutputWorkspace", outWSName));
     TS_ASSERT_THROWS_NOTHING(
@@ -160,7 +161,7 @@ public:
       // Remove workspace from the data service.
       AnalysisDataService::Instance().remove(outWSName);
     }
-    AnalysisDataService::Instance().remove("MDEWS");
+    AnalysisDataService::Instance().remove("MDWS");
   }
 
   /** Running the algo twice with same output workspace = replace the output,
@@ -225,14 +226,14 @@ public:
     FrameworkManager::Instance().exec(
         "BinMD", 14, "AxisAligned", "1", "AlignedDim0", "Q_lab_x,-10,10,100",
         "AlignedDim1", "Q_lab_y,-10,10,100", "AlignedDim2",
-        "Q_lab_z,-10,10,100", "IterateEvents", "1", "InputWorkspace", "MDEWS",
-        "OutputWorkspace", "MDEWS");
+        "Q_lab_z,-10,10,100", "IterateEvents", "1", "InputWorkspace", "MDWS",
+        "OutputWorkspace", "MDWS");
 
     FindPeaksMD alg;
     alg.setRethrows(true);
     TS_ASSERT_THROWS_NOTHING(alg.initialize());
     TS_ASSERT(alg.isInitialized());
-    TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("InputWorkspace", "MDEWS"));
+    TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("InputWorkspace", "MDWS"));
     TS_ASSERT_THROWS_NOTHING(
         alg.setPropertyValue("OutputWorkspace", "place_holder"));
     TS_ASSERT_THROWS_NOTHING(
@@ -245,7 +246,113 @@ public:
     TS_ASSERT_THROWS_NOTHING(alg.setProperty("SignalThresholdFactor", 1.3));
     TS_ASSERT_THROWS(alg.execute(), const std::runtime_error &);
 
-    AnalysisDataService::Instance().remove("MDEWS");
+    AnalysisDataService::Instance().remove("MDWS");
+  }
+
+  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",
+        "MDWS");
+
+    if (expInfo) {
+      Instrument_sptr inst =
+          ComponentCreationHelper::createTestInstrumentRectangular2(1, 100,
+                                                                    0.05);
+      IMDEventWorkspace_sptr ws;
+      TS_ASSERT_THROWS_NOTHING(
+          ws = AnalysisDataService::Instance().retrieveWS<IMDEventWorkspace>(
+              "MDWS"));
+      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", "MDWS", "OutputWorkspace", "MDWS");
+    }
+
+    std::string outWSName("peaksFound");
+    FindPeaksMD alg;
+    TS_ASSERT_THROWS_NOTHING(alg.initialize())
+    TS_ASSERT(alg.isInitialized())
+    TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("InputWorkspace", "MDWS"));
+    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("OutputType", "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("MDWS");
+  }
+
+  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);
   }
 };
 
@@ -286,7 +393,7 @@ public:
     FindPeaksMD alg;
     alg.initialize();
 
-    alg.setPropertyValue("InputWorkspace", "MDEWS");
+    alg.setPropertyValue("InputWorkspace", "MDWS");
     alg.setPropertyValue("OutputWorkspace", outWSName);
     alg.setPropertyValue("DensityThresholdFactor", "2.0");
     alg.setPropertyValue("PeakDistanceThreshold", "0.7");

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,
+                                                     OutputType='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):

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',
+                    OutputType='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')

docs/source/algorithms/FindPeaksMD-v1.rst

@@ -33,9 +33,13 @@ The algorithm proceeds in this way:
 
 -  This is repeated until we find up to MaxPeaks peaks.
 
-Each peak created is placed in the output
-:ref:`PeaksWorkspace <PeaksWorkspace>`, which can be a new workspace or
-replace the old one.
+Each peak created is placed in the output :ref:`PeaksWorkspace
+<PeaksWorkspace>` or :ref:`LeanElasticPeaksWorkspace
+<LeanElasticPeaksWorkspace>` (depending on the `OutputType` option),
+which can be a new workspace or replace the old one. If `OutputType`
+is the default `Automatic` then the output type will be PeakWorkspace
+unless the input workspace doesn't contain an experiment info in which
+case it will default to LeanElasticPeaksWorkspace.
 
 This algorithm works on a :ref:`MDHistoWorkspace <MDHistoWorkspace>`
 resulting from the :ref:`algm-BinMD` algorithm also. It works in the

docs/source/concepts/LeanElasticPeaksWorkspace.rst

@@ -10,6 +10,7 @@ of single crystal LeanElasticPeak objects. It is the equivalent to the
 Creating a LeanElasticPeaksWorkspace
 ------------------------------------
 
+* :ref:`FindPeaksMD <algm-FindPeaksMD>` will find peaks in reciprocal space in a :ref:`MDWorkspace <MDWorkspace>`.
 * :ref:`CreatePeaksWorkspace <algm-CreatePeaksWorkspace>` will create an empty LeanElasticPeaksWorkspace that you can then edit.
 
 Viewing a LeanElasticPeaksWorkspace

docs/source/release/v6.1.0/diffraction.rst

@@ -99,7 +99,8 @@ this peak is a Q-sample vector. There are a number of modifications
 made to facilitate this.
 
 - New LeanElasticPeak and LeanElasticPeakWorkspace has been created :ref:`LeanElasticPeaksWorkspace <LeanElasticPeaksWorkspace>`
-- :ref:`CreatePeaksWorkspace <algm-CreatePeaksWorkspace>` has been modified to optionally create a  :ref:`LeanElasticPeaksWorkspace <LeanElasticPeaksWorkspace>`.
+- :ref:`CreatePeaksWorkspace <algm-CreatePeaksWorkspace>` has been modified to optionally create a :ref:`LeanElasticPeaksWorkspace <LeanElasticPeaksWorkspace>`.
+- :ref:`FindPeaksMD <algm-FindPeaksMD>` has been modified to optionally create a :ref:`LeanElasticPeaksWorkspace <LeanElasticPeaksWorkspace>`.
 - These following other algorithms have either been made to work or confirmed to already work with the LeanElasticPeak:
 
    - :ref:`algm-AddPeakHKL`
@@ -113,6 +114,8 @@ made to facilitate this.
    - :ref:`algm-FindUBUsingMinMaxD`
    - :ref:`algm-IndexPeaks`
    - :ref:`algm-IntegratePeaksMD`
+   - :ref:`algm-LoadNexusProcessed`
+   - :ref:`algm-SaveNexusProcessed`
    - :ref:`algm-SelectCellOfType`
    - :ref:`algm-SelectCellWithForm`
    - :ref:`algm-ShowPossibleCells`