diff --git a/Framework/MDAlgorithms/inc/MantidMDAlgorithms/FindPeaksMD.h b/Framework/MDAlgorithms/inc/MantidMDAlgorithms/FindPeaksMD.h
index f152fb820814381fa219ee40cd7cb4942b241e9c..64dc5bb72c26030d3c3a0a0bb1536cf1990daeb4 100644
--- a/Framework/MDAlgorithms/inc/MantidMDAlgorithms/FindPeaksMD.h
+++ b/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
diff --git a/Framework/MDAlgorithms/src/FindPeaksMD.cpp b/Framework/MDAlgorithms/src/FindPeaksMD.cpp
index 22f5f9ac99dafd87b7b2bdab07c75939e8384014..3b4418be39eb815de36ba48dd95955350a53ceb1 100644
--- a/Framework/MDAlgorithms/src/FindPeaksMD.cpp
+++ b/Framework/MDAlgorithms/src/FindPeaksMD.cpp
@@ -8,7 +8,6 @@
#include "MantidAPI/Run.h"
#include "MantidDataObjects/MDEventFactory.h"
#include "MantidDataObjects/MDHistoWorkspace.h"
-#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidGeometry/Crystal/EdgePixel.h"
#include "MantidGeometry/Instrument/Goniometer.h"
#include "MantidGeometry/Objects/InstrumentRayTracer.h"
@@ -227,8 +226,11 @@ void FindPeaksMD::init() {
std::make_unique<EnabledWhenProperty>(
"CalculateGoniometerForCW",
Mantid::Kernel::ePropertyCriterion::IS_NOT_DEFAULT));
-
- declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace>>(
+ std::vector<std::string> peakTypes = {"Automatic", "Peak", "LeanElasticPeak"};
+ declareProperty("OutputType", "Automatic",
+ std::make_shared<StringListValidator>(peakTypes),
+ "Type of Peak in OutputWorkspace");
+ declareProperty(std::make_unique<WorkspaceProperty<IPeaksWorkspace>>(
"OutputWorkspace", "", Direction::Output),
"An output PeaksWorkspace with the peaks' found positions.");
@@ -245,13 +247,25 @@ void FindPeaksMD::init() {
//----------------------------------------------------------------------------------------------
/** Extract needed data from the workspace's experiment info */
-void FindPeaksMD::readExperimentInfo(const ExperimentInfo_sptr &ei,
- const IMDWorkspace_sptr &ws) {
+void FindPeaksMD::readExperimentInfo(const ExperimentInfo_sptr &ei) {
// Instrument associated with workspace
inst = ei->getInstrument();
// Find the run number
m_runNumber = ei->getRunNumber();
+ // Find the goniometer rotation matrix
+ m_goniometer =
+ Mantid::Kernel::Matrix<double>(3, 3, true); // Default IDENTITY matrix
+ try {
+ m_goniometer = ei->mutableRun().getGoniometerMatrix();
+ } catch (std::exception &e) {
+ g_log.warning() << "Error finding goniometer matrix. It will not be set in "
+ "the peaks found.\n";
+ g_log.warning() << e.what() << '\n';
+ }
+}
+
+void FindPeaksMD::checkWorkspaceDims(const IMDWorkspace_sptr &ws) {
// Check that the workspace dimensions are in Q-sample-frame or Q-lab-frame.
std::string dim0 = ws->getDimension(0)->getName();
if (dim0 == "H") {
@@ -265,16 +279,22 @@ void FindPeaksMD::readExperimentInfo(const ExperimentInfo_sptr &ei,
else
throw std::runtime_error(
"Unexpected dimensions: need either Q_lab_x or Q_sample_x.");
+}
- // Find the goniometer rotation matrix
- m_goniometer =
- Mantid::Kernel::Matrix<double>(3, 3, true); // Default IDENTITY matrix
- try {
- m_goniometer = ei->mutableRun().getGoniometerMatrix();
- } catch (std::exception &e) {
- g_log.warning() << "Error finding goniometer matrix. It will not be set in "
- "the peaks found.\n";
- g_log.warning() << e.what() << '\n';
+void FindPeaksMD::determineOutputType(const std::string peakType,
+ const uint16_t numExperimentInfo) {
+ // This method will be expanded later to check a property on the
+ // input workspace which can specify a default peak type for that
+ // instrument.
+ m_leanElasticPeak = false;
+ if (peakType == "Automatic") {
+ if (numExperimentInfo == 0)
+ m_leanElasticPeak = true;
+ } else if (peakType == "LeanElasticPeak") {
+ m_leanElasticPeak = true;
+ } else { // Peak
+ if (numExperimentInfo == 0)
+ throw std::runtime_error("Cannot create Peak output with 0 expInfo");
}
}
@@ -301,6 +321,20 @@ void FindPeaksMD::addPeak(const V3D &Q, const double binCount,
}
}
+//----------------------------------------------------------------------------------------------
+/** Create and add a LeanElasticPeak to the output workspace
+ *
+ * @param Q :: Q_lab or Q_sample, depending on workspace
+ * @param binCount :: bin count to give to the peak.
+ * @param useGoniometer :: if to include set goniometer from the input workspace
+ * in the peak
+ */
+void FindPeaksMD::addLeanElasticPeak(const V3D &Q, const double binCount,
+ const bool useGoniometer) {
+ auto p = this->createLeanElasticPeak(Q, binCount, useGoniometer);
+ peakWS->addPeak(*p);
+}
+
/**
* Creates a Peak object from Q & bin count
* */
@@ -324,9 +358,9 @@ FindPeaksMD::createPeak(const Mantid::Kernel::V3D &Q, const double binCount,
if (inst->hasParameter("wavelength")) {
wavelength = inst->getNumberParameter("wavelength").at(0);
} else {
- throw std::runtime_error(
- "Could not get wavelength, neither Wavelength algorithm property "
- "set nor instrument wavelength parameter");
+ throw std::runtime_error("Could not get wavelength, neither "
+ "Wavelength algorithm property "
+ "set nor instrument wavelength parameter");
}
}
@@ -359,9 +393,32 @@ FindPeaksMD::createPeak(const Mantid::Kernel::V3D &Q, const double binCount,
return p;
}
+/**
+ * Creates a Peak object from Q & bin count
+ * */
+std::shared_ptr<DataObjects::LeanElasticPeak>
+FindPeaksMD::createLeanElasticPeak(const Mantid::Kernel::V3D &Q,
+ const double binCount,
+ const bool useGoniometer) {
+ std::shared_ptr<DataObjects::LeanElasticPeak> p;
+ if (dimType == QSAMPLE) {
+ if (useGoniometer)
+ p = std::make_shared<LeanElasticPeak>(Q, m_goniometer);
+ else
+ p = std::make_shared<LeanElasticPeak>(Q);
+ } else {
+ throw std::invalid_argument(
+ "Cannot find peaks unless the dimension is QSAMPLE");
+ }
+
+ p->setBinCount(binCount);
+ // Save the run number found before.
+ p->setRunNumber(m_runNumber);
+ return p;
+}
//----------------------------------------------------------------------------------------------
-/** Integrate the peaks of the workspace using parameters saved in the algorithm
- * class
+/** Integrate the peaks of the workspace using parameters saved in the
+ * algorithm class
* @param ws :: MDEventWorkspace to integrate
*/
template <typename MDE, size_t nd>
@@ -384,168 +441,127 @@ void FindPeaksMD::findPeaks(typename MDEventWorkspace<MDE, nd>::sptr ws) {
// Make sure all centroids are fresh
// ws->getBox()->refreshCentroid();
- uint16_t nexp = ws->getNumExperimentInfo();
-
- if (nexp == 0)
- throw std::runtime_error(
- "No instrument was found in the MDEventWorkspace. Cannot find peaks.");
-
- for (uint16_t iexp = 0; iexp < ws->getNumExperimentInfo(); iexp++) {
- ExperimentInfo_sptr ei = ws->getExperimentInfo(iexp);
- this->readExperimentInfo(ei, std::dynamic_pointer_cast<IMDWorkspace>(ws));
-
- Geometry::InstrumentRayTracer tracer(inst);
- // Copy the instrument, sample, run to the peaks workspace.
- peakWS->copyExperimentInfoFrom(ei.get());
-
- // Calculate a threshold below which a box is too diffuse to be considered a
- // peak.
- signal_t threshold =
- m_useNumberOfEventsNormalization
- ? ws->getBox()->getSignalByNEvents() * m_signalThresholdFactor
- : ws->getBox()->getSignalNormalized() * DensityThresholdFactor;
- threshold *= m_densityScaleFactor;
-
- if (!std::isfinite(threshold)) {
- g_log.warning()
- << "Infinite or NaN overall density found. Your input data "
- "may be invalid. Using a 0 threshold instead.\n";
- threshold = 0;
- }
- g_log.information() << "Threshold signal density: " << threshold << '\n';
-
- using boxPtr = API::IMDNode *;
- // We will fill this vector with pointers to all the boxes (up to a given
- // depth)
- typename std::vector<API::IMDNode *> boxes;
-
- // Get all the MDboxes
- progress(0.10, "Getting Boxes");
- ws->getBox()->getBoxes(boxes, 1000, true);
-
- // This pair is the <density, ptr to the box>
- using dens_box = std::pair<double, API::IMDNode *>;
-
- // Map that will sort the boxes by increasing density. The key = density;
- // value = box *.
- typename std::multimap<double, API::IMDNode *> sortedBoxes;
-
- // --------------- Sort and Filter by Density -----------------------------
- progress(0.20, "Sorting Boxes by Density");
- auto it1 = boxes.begin();
- auto it1_end = boxes.end();
- for (; it1 != it1_end; it1++) {
- auto box = *it1;
- double value = m_useNumberOfEventsNormalization
- ? box->getSignalByNEvents()
- : box->getSignalNormalized();
- value *= m_densityScaleFactor;
- // Skip any boxes with too small a signal value.
- if (value > threshold)
- sortedBoxes.insert(dens_box(value, box));
- }
+ // Calculate a threshold below which a box is too diffuse to be considered a
+ // peak.
+ signal_t threshold =
+ m_useNumberOfEventsNormalization
+ ? ws->getBox()->getSignalByNEvents() * m_signalThresholdFactor
+ : ws->getBox()->getSignalNormalized() * DensityThresholdFactor;
+ threshold *= m_densityScaleFactor;
+
+ if (!std::isfinite(threshold)) {
+ g_log.warning() << "Infinite or NaN overall density found. Your input data "
+ "may be invalid. Using a 0 threshold instead.\n";
+ threshold = 0;
+ }
+ g_log.information() << "Threshold signal density: " << threshold << '\n';
+
+ using boxPtr = API::IMDNode *;
+ // We will fill this vector with pointers to all the boxes (up to a given
+ // depth)
+ typename std::vector<API::IMDNode *> boxes;
+
+ // Get all the MDboxes
+ progress(0.10, "Getting Boxes");
+ ws->getBox()->getBoxes(boxes, 1000, true);
+
+ // This pair is the <density, ptr to the box>
+ using dens_box = std::pair<double, API::IMDNode *>;
+
+ // Map that will sort the boxes by increasing density. The key = density;
+ // value = box *.
+ typename std::multimap<double, API::IMDNode *> sortedBoxes;
+
+ // --------------- Sort and Filter by Density -----------------------------
+ progress(0.20, "Sorting Boxes by Density");
+ auto it1 = boxes.begin();
+ auto it1_end = boxes.end();
+ for (; it1 != it1_end; it1++) {
+ auto box = *it1;
+ double value = m_useNumberOfEventsNormalization
+ ? box->getSignalByNEvents()
+ : box->getSignalNormalized();
+ value *= m_densityScaleFactor;
+ // Skip any boxes with too small a signal value.
+ if (value > threshold)
+ sortedBoxes.insert(dens_box(value, box));
+ }
- // --------------- Find Peak Boxes -----------------------------
- // List of chosen possible peak boxes.
- std::vector<API::IMDNode *> peakBoxes;
+ // --------------- Find Peak Boxes -----------------------------
+ // List of chosen possible peak boxes.
+ std::vector<API::IMDNode *> peakBoxes;
- prog = std::make_unique<Progress>(this, 0.30, 0.95, m_maxPeaks);
+ prog = std::make_unique<Progress>(this, 0.30, 0.95, m_maxPeaks);
- // used for selecting method for calculating BinCount
- bool isMDEvent(ws->id().find("MDEventWorkspace") != std::string::npos);
+ // used for selecting method for calculating BinCount
+ bool isMDEvent(ws->id().find("MDEventWorkspace") != std::string::npos);
- int64_t numBoxesFound = 0;
- // Now we go (backwards) through the map
- // e.g. from highest density down to lowest density.
- typename std::multimap<double, boxPtr>::reverse_iterator it2;
- auto it2_end = sortedBoxes.rend();
- for (it2 = sortedBoxes.rbegin(); it2 != it2_end; ++it2) {
- signal_t density = it2->first;
- boxPtr box = it2->second;
+ int64_t numBoxesFound = 0;
+ // Now we go (backwards) through the map
+ // e.g. from highest density down to lowest density.
+ typename std::multimap<double, boxPtr>::reverse_iterator it2;
+ auto it2_end = sortedBoxes.rend();
+ for (it2 = sortedBoxes.rbegin(); it2 != it2_end; ++it2) {
+ signal_t density = it2->first;
+ boxPtr box = it2->second;
#ifndef MDBOX_TRACK_CENTROID
- coord_t boxCenter[nd];
- box->calculateCentroid(boxCenter);
+ coord_t boxCenter[nd];
+ box->calculateCentroid(boxCenter);
#else
- const coord_t *boxCenter = box->getCentroid();
+ const coord_t *boxCenter = box->getCentroid();
#endif
- // Compare to all boxes already picked.
- bool badBox = false;
- for (auto &peakBoxe : peakBoxes) {
+ // Compare to all boxes already picked.
+ bool badBox = false;
+ for (auto &peakBoxe : peakBoxes) {
#ifndef MDBOX_TRACK_CENTROID
- coord_t otherCenter[nd];
- (*it3)->calculateCentroid(otherCenter);
+ coord_t otherCenter[nd];
+ (*it3)->calculateCentroid(otherCenter);
#else
- const coord_t *otherCenter = peakBoxe->getCentroid();
+ const coord_t *otherCenter = peakBoxe->getCentroid();
#endif
- // Distance between this box and a box we already put in.
- coord_t distSquared = 0.0;
- for (size_t d = 0; d < nd; d++) {
- coord_t dist = otherCenter[d] - boxCenter[d];
- distSquared += (dist * dist);
- }
-
- // Reject this box if it is too close to another previously found box.
- if (distSquared < peakRadiusSquared) {
- badBox = true;
- break;
- }
+ // Distance between this box and a box we already put in.
+ coord_t distSquared = 0.0;
+ for (size_t d = 0; d < nd; d++) {
+ coord_t dist = otherCenter[d] - boxCenter[d];
+ distSquared += (dist * dist);
}
- // The box was not rejected for another reason.
- if (!badBox) {
- if (numBoxesFound++ >= m_maxPeaks) {
- g_log.notice() << "Number of peaks found exceeded the limit of "
- << m_maxPeaks << ". Stopping peak finding.\n";
- break;
- }
+ // Reject this box if it is too close to another previously found box.
+ if (distSquared < peakRadiusSquared) {
+ badBox = true;
+ break;
+ }
+ }
- peakBoxes.emplace_back(box);
- g_log.debug() << "Found box at ";
- for (size_t d = 0; d < nd; d++)
- g_log.debug() << (d > 0 ? "," : "") << boxCenter[d];
- g_log.debug() << "; Density = " << density << '\n';
- // Report progres for each box found.
- prog->report("Finding Peaks");
+ // The box was not rejected for another reason.
+ if (!badBox) {
+ if (numBoxesFound++ >= m_maxPeaks) {
+ g_log.notice() << "Number of peaks found exceeded the limit of "
+ << m_maxPeaks << ". Stopping peak finding.\n";
+ break;
}
+
+ peakBoxes.emplace_back(box);
+ g_log.debug() << "Found box at ";
+ for (size_t d = 0; d < nd; d++)
+ g_log.debug() << (d > 0 ? "," : "") << boxCenter[d];
+ g_log.debug() << "; Density = " << density << '\n';
+ // Report progres for each box found.
+ prog->report("Finding Peaks");
}
+ }
- prog->resetNumSteps(numBoxesFound, 0.95, 1.0);
+ prog->resetNumSteps(numBoxesFound, 0.95, 1.0);
+ uint16_t numExperimentInfo = ws->getNumExperimentInfo();
+
+ if (numExperimentInfo == 0) {
// --- Convert the "boxes" to peaks ----
for (auto box : peakBoxes) {
- // If no events from this experimental contribute to the box then skip
- if (nexp > 1) {
- auto *mdbox = dynamic_cast<MDBox<MDE, nd> *>(box);
- const std::vector<MDE> &events = mdbox->getEvents();
- if (std::none_of(events.cbegin(), events.cend(),
- [&iexp, &nexp](MDE event) {
- return event.getRunIndex() == iexp ||
- event.getRunIndex() >= nexp;
- }))
- continue;
- }
-
- // If multiple goniometers than use the average one from the
- // events in the box, that matches this runIndex, this assumes
- // the events are added in same order as the goniometers
- if (ei->run().getNumGoniometers() > 1) {
- const std::vector<MDE> &events =
- dynamic_cast<MDBox<MDE, nd> *>(box)->getEvents();
- double sum = 0;
- double count = 0;
- for (const auto &event : events) {
- if (event.getRunIndex() == iexp) {
- sum += event.getGoniometerIndex();
- count++;
- }
- }
- m_goniometer = ei->mutableRun().getGoniometerMatrix(lrint(sum / count));
- }
// The center of the box = Q in the lab frame
-
#ifndef MDBOX_TRACK_CENTROID
coord_t boxCenter[nd];
box->calculateCentroid(boxCenter);
@@ -556,42 +572,111 @@ void FindPeaksMD::findPeaks(typename MDEventWorkspace<MDE, nd>::sptr ws) {
// Q of the centroid of the box
V3D Q(boxCenter[0], boxCenter[1], boxCenter[2]);
- // The "bin count" used will be the box density or the number of events in
- // the box
+ // The "bin count" used will be the box density.
double binCount = box->getSignalNormalized() * m_densityScaleFactor;
if (isMDEvent)
binCount = static_cast<double>(box->getNPoints());
- try {
- auto p = this->createPeak(Q, binCount, tracer);
- if (m_addDetectors) {
- auto mdBox = dynamic_cast<MDBoxBase<MDE, nd> *>(box);
- if (!mdBox) {
- throw std::runtime_error("Failed to cast box to MDBoxBase");
+ // Create the peak
+ addLeanElasticPeak(Q, binCount);
+
+ // Report progres for each box found.
+ prog->report("Adding Peaks");
+
+ } // for each box found
+ } else {
+ for (uint16_t iexp = 0; iexp < ws->getNumExperimentInfo(); iexp++) {
+ ExperimentInfo_sptr ei = ws->getExperimentInfo(iexp);
+ this->readExperimentInfo(ei);
+
+ Geometry::InstrumentRayTracer tracer(inst);
+ // Copy the instrument, sample, run to the peaks workspace.
+ peakWS->copyExperimentInfoFrom(ei.get());
+
+ // --- Convert the "boxes" to peaks ----
+ for (auto box : peakBoxes) {
+ // If no events from this experimental contribute to the box then skip
+ if (numExperimentInfo > 1) {
+ auto *mdbox = dynamic_cast<MDBox<MDE, nd> *>(box);
+ const std::vector<MDE> &events = mdbox->getEvents();
+ if (std::none_of(events.cbegin(), events.cend(),
+ [&iexp, &numExperimentInfo](MDE event) {
+ return event.getRunIndex() == iexp ||
+ event.getRunIndex() >= numExperimentInfo;
+ }))
+ continue;
+ }
+
+ // If multiple goniometers than use the average one from the
+ // events in the box, that matches this runIndex, this assumes
+ // the events are added in same order as the goniometers
+ if (ei->run().getNumGoniometers() > 1) {
+ const std::vector<MDE> &events =
+ dynamic_cast<MDBox<MDE, nd> *>(box)->getEvents();
+ double sum = 0;
+ double count = 0;
+ for (const auto &event : events) {
+ if (event.getRunIndex() == iexp) {
+ sum += event.getGoniometerIndex();
+ count++;
+ }
}
- addDetectors(*p, *mdBox);
+ m_goniometer =
+ ei->mutableRun().getGoniometerMatrix(lrint(sum / count));
}
- if (p->getDetectorID() != -1) {
- if (m_edge > 0) {
- if (!edgePixel(inst, p->getBankName(), p->getCol(), p->getRow(),
- m_edge))
- peakWS->addPeak(*p);
- ;
- } else {
- peakWS->addPeak(*p);
+ // The center of the box = Q in the lab frame
+
+#ifndef MDBOX_TRACK_CENTROID
+ coord_t boxCenter[nd];
+ box->calculateCentroid(boxCenter);
+#else
+ const coord_t *boxCenter = box->getCentroid();
+#endif
+
+ // Q of the centroid of the box
+ V3D Q(boxCenter[0], boxCenter[1], boxCenter[2]);
+
+ // The "bin count" used will be the box density or the number of events
+ // in the box
+ double binCount = box->getSignalNormalized() * m_densityScaleFactor;
+ if (isMDEvent)
+ binCount = static_cast<double>(box->getNPoints());
+
+ if (m_leanElasticPeak) {
+ addLeanElasticPeak(Q, binCount, true);
+ } else {
+ try {
+ auto p = this->createPeak(Q, binCount, tracer);
+ if (m_addDetectors) {
+ auto mdBox = dynamic_cast<MDBoxBase<MDE, nd> *>(box);
+ if (!mdBox) {
+ throw std::runtime_error("Failed to cast box to MDBoxBase");
+ }
+ addDetectors(*p, *mdBox);
+ }
+ if (p->getDetectorID() != -1) {
+ if (m_edge > 0) {
+ if (!edgePixel(inst, p->getBankName(), p->getCol(), p->getRow(),
+ m_edge))
+ peakWS->addPeak(*p);
+ ;
+ } else {
+ peakWS->addPeak(*p);
+ }
+ g_log.information() << "Add new peak with Q-center = " << Q[0]
+ << ", " << Q[1] << ", " << Q[2] << "\n";
+ }
+ } catch (std::exception &e) {
+ g_log.notice() << "Error creating peak at " << Q << " because of '"
+ << e.what() << "'. Peak will be skipped.\n";
}
- g_log.information() << "Add new peak with Q-center = " << Q[0] << ", "
- << Q[1] << ", " << Q[2] << "\n";
}
- } catch (std::exception &e) {
- g_log.notice() << "Error creating peak at " << Q << " because of '"
- << e.what() << "'. Peak will be skipped.\n";
- }
- // Report progress for each box found.
- prog->report("Adding Peaks");
+ // Report progress for each box found.
+ prog->report("Adding Peaks");
- } // for each box found
+ } // for each box found
+ }
}
g_log.notice() << "Number of peaks found: " << peakWS->getNumberPeaks()
<< '\n';
@@ -611,105 +696,94 @@ void FindPeaksMD::findPeaksHisto(
g_log.warning("Workspace is an MDHistoWorkspace. Resultant PeaksWorkspaces "
"will not contain full detector information.");
- if (ws->getNumExperimentInfo() == 0)
- throw std::runtime_error(
- "No instrument was found in the workspace. Cannot find peaks.");
-
- for (uint16_t iexp = 0; iexp < ws->getNumExperimentInfo(); iexp++) {
- ExperimentInfo_sptr ei = ws->getExperimentInfo(iexp);
- this->readExperimentInfo(ei, std::dynamic_pointer_cast<IMDWorkspace>(ws));
- Geometry::InstrumentRayTracer tracer(inst);
-
- // Copy the instrument, sample, run to the peaks workspace.
- peakWS->copyExperimentInfoFrom(ei.get());
-
- // This pair is the <density, box index>
- using dens_box = std::pair<double, size_t>;
-
- // Map that will sort the boxes by increasing density. The key = density;
- // value = box index.
- std::multimap<double, size_t> sortedBoxes;
-
- size_t numBoxes = ws->getNPoints();
-
- // --------- Count the overall signal density -----------------------------
- progress(0.10, "Counting Total Signal");
- double totalSignal = 0;
- for (size_t i = 0; i < numBoxes; i++)
- totalSignal += ws->getSignalAt(i);
- // Calculate the threshold density
- double thresholdDensity =
- (totalSignal * ws->getInverseVolume() / double(numBoxes)) *
- DensityThresholdFactor * m_densityScaleFactor;
- if (!std::isfinite(thresholdDensity)) {
- g_log.warning()
- << "Infinite or NaN overall density found. Your input data "
- "may be invalid. Using a 0 threshold instead.\n";
- thresholdDensity = 0;
- }
- g_log.information() << "Threshold signal density: " << thresholdDensity
- << '\n';
-
- // -------------- Sort and Filter by Density -----------------------------
- progress(0.20, "Sorting Boxes by Density");
- for (size_t i = 0; i < numBoxes; i++) {
- double density = ws->getSignalNormalizedAt(i) * m_densityScaleFactor;
- // Skip any boxes with too small a signal density.
- if (density > thresholdDensity)
- sortedBoxes.insert(dens_box(density, i));
- }
-
- // --------------- Find Peak Boxes -----------------------------
- // List of chosen possible peak boxes.
- std::vector<size_t> peakBoxes;
-
- prog = std::make_unique<Progress>(this, 0.30, 0.95, m_maxPeaks);
-
- int64_t numBoxesFound = 0;
- // Now we go (backwards) through the map
- // e.g. from highest density down to lowest density.
- std::multimap<double, size_t>::reverse_iterator it2;
- auto it2_end = sortedBoxes.rend();
- for (it2 = sortedBoxes.rbegin(); it2 != it2_end; ++it2) {
- signal_t density = it2->first;
- size_t index = it2->second;
- // Get the center of the box
- VMD boxCenter = ws->getCenter(index);
-
- // Compare to all boxes already picked.
- bool badBox = false;
- for (auto &peakBoxe : peakBoxes) {
- VMD otherCenter = ws->getCenter(peakBoxe);
-
- // Distance between this box and a box we already put in.
- coord_t distSquared = 0.0;
- for (size_t d = 0; d < nd; d++) {
- coord_t dist = otherCenter[d] - boxCenter[d];
- distSquared += (dist * dist);
- }
+ // This pair is the <density, box index>
+ using dens_box = std::pair<double, size_t>;
+
+ // Map that will sort the boxes by increasing density. The key = density;
+ // value = box index.
+ std::multimap<double, size_t> sortedBoxes;
+
+ size_t numBoxes = ws->getNPoints();
+
+ // --------- Count the overall signal density -----------------------------
+ progress(0.10, "Counting Total Signal");
+ double totalSignal = 0;
+ for (size_t i = 0; i < numBoxes; i++)
+ totalSignal += ws->getSignalAt(i);
+ // Calculate the threshold density
+ double thresholdDensity =
+ (totalSignal * ws->getInverseVolume() / double(numBoxes)) *
+ DensityThresholdFactor * m_densityScaleFactor;
+ if (!std::isfinite(thresholdDensity)) {
+ g_log.warning() << "Infinite or NaN overall density found. Your input data "
+ "may be invalid. Using a 0 threshold instead.\n";
+ thresholdDensity = 0;
+ }
+ g_log.information() << "Threshold signal density: " << thresholdDensity
+ << '\n';
+
+ // -------------- Sort and Filter by Density -----------------------------
+ progress(0.20, "Sorting Boxes by Density");
+ for (size_t i = 0; i < numBoxes; i++) {
+ double density = ws->getSignalNormalizedAt(i) * m_densityScaleFactor;
+ // Skip any boxes with too small a signal density.
+ if (density > thresholdDensity)
+ sortedBoxes.insert(dens_box(density, i));
+ }
- // Reject this box if it is too close to another previously found box.
- if (distSquared < peakRadiusSquared) {
- badBox = true;
- break;
- }
+ // --------------- Find Peak Boxes -----------------------------
+ // List of chosen possible peak boxes.
+ std::vector<size_t> peakBoxes;
+
+ prog = std::make_unique<Progress>(this, 0.30, 0.95, m_maxPeaks);
+
+ int64_t numBoxesFound = 0;
+ // Now we go (backwards) through the map
+ // e.g. from highest density down to lowest density.
+ std::multimap<double, size_t>::reverse_iterator it2;
+ auto it2_end = sortedBoxes.rend();
+ for (it2 = sortedBoxes.rbegin(); it2 != it2_end; ++it2) {
+ signal_t density = it2->first;
+ size_t index = it2->second;
+ // Get the center of the box
+ VMD boxCenter = ws->getCenter(index);
+
+ // Compare to all boxes already picked.
+ bool badBox = false;
+ for (auto &peakBoxe : peakBoxes) {
+ VMD otherCenter = ws->getCenter(peakBoxe);
+
+ // Distance between this box and a box we already put in.
+ coord_t distSquared = 0.0;
+ for (size_t d = 0; d < nd; d++) {
+ coord_t dist = otherCenter[d] - boxCenter[d];
+ distSquared += (dist * dist);
}
- // The box was not rejected for another reason.
- if (!badBox) {
- if (numBoxesFound++ >= m_maxPeaks) {
- g_log.notice() << "Number of peaks found exceeded the limit of "
- << m_maxPeaks << ". Stopping peak finding.\n";
- break;
- }
+ // Reject this box if it is too close to another previously found box.
+ if (distSquared < peakRadiusSquared) {
+ badBox = true;
+ break;
+ }
+ }
- peakBoxes.emplace_back(index);
- g_log.debug() << "Found box at index " << index;
- g_log.debug() << "; Density = " << density << '\n';
- // Report progres for each box found.
- prog->report("Finding Peaks");
+ // The box was not rejected for another reason.
+ if (!badBox) {
+ if (numBoxesFound++ >= m_maxPeaks) {
+ g_log.notice() << "Number of peaks found exceeded the limit of "
+ << m_maxPeaks << ". Stopping peak finding.\n";
+ break;
}
+
+ peakBoxes.emplace_back(index);
+ g_log.debug() << "Found box at index " << index;
+ g_log.debug() << "; Density = " << density << '\n';
+ // Report progres for each box found.
+ prog->report("Finding Peaks");
}
+ }
+
+ if (ws->getNumExperimentInfo() == 0) {
// --- Convert the "boxes" to peaks ----
for (auto index : peakBoxes) {
// The center of the box = Q in the lab frame
@@ -722,36 +796,84 @@ void FindPeaksMD::findPeaksHisto(
double binCount = ws->getSignalNormalizedAt(index) * m_densityScaleFactor;
// Create the peak
- addPeak(Q, binCount, tracer);
+ addLeanElasticPeak(Q, binCount);
// Report progres for each box found.
prog->report("Adding Peaks");
} // for each box found
+ } else {
+ for (uint16_t iexp = 0; iexp < ws->getNumExperimentInfo(); iexp++) {
+ ExperimentInfo_sptr ei = ws->getExperimentInfo(iexp);
+ this->readExperimentInfo(ei);
+ Geometry::InstrumentRayTracer tracer(inst);
+
+ // Copy the instrument, sample, run to the peaks workspace.
+ peakWS->copyExperimentInfoFrom(ei.get());
+
+ // --- Convert the "boxes" to peaks ----
+ for (auto index : peakBoxes) {
+ // The center of the box = Q in the lab frame
+ VMD boxCenter = ws->getCenter(index);
+
+ // Q of the centroid of the box
+ V3D Q(boxCenter[0], boxCenter[1], boxCenter[2]);
+
+ // The "bin count" used will be the box density.
+ double binCount =
+ ws->getSignalNormalizedAt(index) * m_densityScaleFactor;
+
+ // Create the peak
+ if (m_leanElasticPeak)
+ addLeanElasticPeak(Q, binCount, true);
+ else
+ addPeak(Q, binCount, tracer);
+
+ // Report progres for each box found.
+ prog->report("Adding Peaks");
+
+ } // for each box found
+ }
}
g_log.notice() << "Number of peaks found: " << peakWS->getNumberPeaks()
<< '\n';
-}
+} // namespace MDAlgorithms
//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
*/
void FindPeaksMD::exec() {
- bool AppendPeaks = getProperty("AppendPeaks");
-
- // Output peaks workspace, create if needed
- peakWS = getProperty("OutputWorkspace");
- if (!peakWS || !AppendPeaks)
- peakWS = PeaksWorkspace_sptr(new PeaksWorkspace());
-
// The MDEventWorkspace as input
IMDWorkspace_sptr inWS = getProperty("InputWorkspace");
+ checkWorkspaceDims(inWS);
MDHistoWorkspace_sptr inMDHW =
std::dynamic_pointer_cast<MDHistoWorkspace>(inWS);
IMDEventWorkspace_sptr inMDEW =
std::dynamic_pointer_cast<IMDEventWorkspace>(inWS);
+ bool AppendPeaks = getProperty("AppendPeaks");
+
+ uint16_t numExperimentInfo = 0;
+ if (inMDHW)
+ numExperimentInfo = inMDHW->getNumExperimentInfo();
+ else if (inMDEW)
+ numExperimentInfo = inMDEW->getNumExperimentInfo();
+
+ std::string peakType = getProperty("OutputType");
+
+ determineOutputType(peakType, numExperimentInfo);
+
+ // Output peaks workspace, create if needed
+ peakWS = getProperty("OutputWorkspace");
+
+ if (!peakWS || !AppendPeaks) {
+ if (m_leanElasticPeak)
+ peakWS = LeanElasticPeaksWorkspace_sptr(new LeanElasticPeaksWorkspace());
+ else
+ peakWS = PeaksWorkspace_sptr(new PeaksWorkspace());
+ }
+
// Other parameters
double PeakDistanceThreshold = getProperty("PeakDistanceThreshold");
peakRadiusSquared =
@@ -785,7 +907,7 @@ void FindPeaksMD::exec() {
peakWS->sort(criteria);
for (auto i = 0; i != peakWS->getNumberPeaks(); ++i) {
- Peak &p = peakWS->getPeak(i);
+ Mantid::Geometry::IPeak &p = peakWS->getPeak(i);
p.setPeakNumber(i + 1);
}
// Save the output
@@ -805,6 +927,8 @@ std::map<std::string, std::string> FindPeaksMD::validateInputs() {
IMDWorkspace_sptr inWS = getProperty("InputWorkspace");
IMDEventWorkspace_sptr inMDEW =
std::dynamic_pointer_cast<IMDEventWorkspace>(inWS);
+ MDHistoWorkspace_sptr inMDHW =
+ std::dynamic_pointer_cast<MDHistoWorkspace>(inWS);
if (useNumberOfEventsNormalization && !inMDEW) {
result["PeakFindingStrategy"] = "The NumberOfEventsNormalization selection "
@@ -812,6 +936,19 @@ std::map<std::string, std::string> FindPeaksMD::validateInputs() {
"as the input.";
}
+ uint16_t numExperimentInfo = 0;
+ if (inMDHW)
+ numExperimentInfo = inMDHW->getNumExperimentInfo();
+ else if (inMDEW)
+ numExperimentInfo = inMDEW->getNumExperimentInfo();
+
+ std::string peakType = getProperty("OutputType");
+
+ if (peakType == "Peak" && numExperimentInfo == 0)
+ result["OutputType"] =
+ "The InputWorkspace doesn't contain any experiment information so the "
+ "OutputType cannot be Peak.";
+
return result;
}
diff --git a/Framework/MDAlgorithms/test/FindPeaksMDTest.h b/Framework/MDAlgorithms/test/FindPeaksMDTest.h
index 0ba29c344a1658d51272cab7e7d592002e0d4870..a92363972e578950db03b74ce4f56ea4c7586e18 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"
@@ -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");
diff --git a/Testing/SystemTests/tests/framework/ConvertHFIRSCDtoMDETest.py b/Testing/SystemTests/tests/framework/ConvertHFIRSCDtoMDETest.py
index d1bc1ab2b9c77983da9a83d474b89b5b8deb7628..3c2bc4deba30f6c3b7547839ecf01c955f508697 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,
+ 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):
diff --git a/Testing/SystemTests/tests/framework/TOPAZPeakFinding.py b/Testing/SystemTests/tests/framework/TOPAZPeakFinding.py
index 55e6b9354bfd4e1a8921c25df0337e7b5486a8bc..104f80fc754bbae89dee2880ed5317b016ee7bd9 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',
+ 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')
diff --git a/docs/source/algorithms/FindPeaksMD-v1.rst b/docs/source/algorithms/FindPeaksMD-v1.rst
index 2834f2473e5c9df77d6c1fff7f854033d75be3a9..b1a132ae9e850c60566e00977a718db95f9c0a75 100644
--- a/docs/source/algorithms/FindPeaksMD-v1.rst
+++ b/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
diff --git a/docs/source/concepts/LeanElasticPeaksWorkspace.rst b/docs/source/concepts/LeanElasticPeaksWorkspace.rst
index 70a8894e138236d949e92b8152354eb5f3948e86..ca6861724acb369a7736fd8020ed521621d8d3eb 100644
--- a/docs/source/concepts/LeanElasticPeaksWorkspace.rst
+++ b/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
diff --git a/docs/source/release/v6.1.0/diffraction.rst b/docs/source/release/v6.1.0/diffraction.rst
index 3d8a6cdd9545ded42227054b8e0461a3893978eb..b2c5ad09d7754c069a717142a59b40e04a15e862 100644
--- a/docs/source/release/v6.1.0/diffraction.rst
+++ b/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`