From 07fed96619fb1ca758452ebdbe4f0442bfee274d Mon Sep 17 00:00:00 2001
From: Janik Zikovsky <zikovskyjl@ornl.gov>
Date: Wed, 24 Aug 2011 14:25:42 +0000
Subject: [PATCH] Refs #3590: MDEWFindPeaks: AppendPeaks option, replace output
if not checked. Also sorting output by bankname and intensity.
---
.../DataObjects/inc/MantidDataObjects/Peak.h | 2 +
.../inc/MantidDataObjects/PeaksWorkspace.h | 4 +
.../Mantid/Framework/DataObjects/src/Peak.cpp | 47 ++++++++++++
.../Framework/DataObjects/src/PeakColumn.cpp | 26 +------
.../DataObjects/src/PeaksWorkspace.cpp | 74 ++++++++++++++++++-
.../Framework/DataObjects/test/PeakTest.h | 15 ++++
.../DataObjects/test/PeaksWorkspaceTest.h | 46 ++++++++++++
.../Framework/MDEvents/src/MDEWFindPeaks.cpp | 39 +++++++---
.../MDEvents/test/MDEWFindPeaksTest.h | 61 +++++++++++----
9 files changed, 264 insertions(+), 50 deletions(-)
diff --git a/Code/Mantid/Framework/DataObjects/inc/MantidDataObjects/Peak.h b/Code/Mantid/Framework/DataObjects/inc/MantidDataObjects/Peak.h
index d346c62e88f..2e2d7401109 100644
--- a/Code/Mantid/Framework/DataObjects/inc/MantidDataObjects/Peak.h
+++ b/Code/Mantid/Framework/DataObjects/inc/MantidDataObjects/Peak.h
@@ -94,6 +94,8 @@ namespace DataObjects
double getL1() const;
double getL2() const;
+ double getValueByColName(const std::string & name) const;
+
protected:
/// Shared pointer to the instrument (for calculating some values )
Geometry::Instrument_const_sptr m_inst;
diff --git a/Code/Mantid/Framework/DataObjects/inc/MantidDataObjects/PeaksWorkspace.h b/Code/Mantid/Framework/DataObjects/inc/MantidDataObjects/PeaksWorkspace.h
index 1f22e682348..f964e569960 100644
--- a/Code/Mantid/Framework/DataObjects/inc/MantidDataObjects/PeaksWorkspace.h
+++ b/Code/Mantid/Framework/DataObjects/inc/MantidDataObjects/PeaksWorkspace.h
@@ -19,6 +19,8 @@
#include "MantidKernel/System.h"
#include "MantidKernel/V3D.h"
#include <string>
+#include <utility>
+#include <vector>
//IsamplePosition should be IsampleOrientation
@@ -77,6 +79,8 @@ namespace DataObjects
void appendFile( std::string filename, Geometry::Instrument_sptr inst);
+ void sort(std::vector< std::pair<std::string, bool> > & criteria);
+
//---------------------------------------------------------------------------------------------
/** @return the number of peaks
*/
diff --git a/Code/Mantid/Framework/DataObjects/src/Peak.cpp b/Code/Mantid/Framework/DataObjects/src/Peak.cpp
index dde9f75c8a9..425803e28db 100644
--- a/Code/Mantid/Framework/DataObjects/src/Peak.cpp
+++ b/Code/Mantid/Framework/DataObjects/src/Peak.cpp
@@ -2,6 +2,9 @@
#include "MantidGeometry/Instrument/RectangularDetector.h"
#include "MantidKernel/System.h"
#include "MantidGeometry/Objects/InstrumentRayTracer.h"
+#include <algorithm>
+#include <string>
+#include <cctype>
using namespace Mantid;
using namespace Mantid::Kernel;
@@ -655,6 +658,50 @@ namespace DataObjects
return (detPos - samplePos).norm();
}
+ // -------------------------------------------------------------------------------------
+ /** Helper function for displaying/sorting peaks in MantidPlot
+ *
+ * @param name :: name of the column in the table workspace
+ * @return a double representing that value (if that's possible)
+ * @throw std::runtime_error if you asked for a column that can't convert to double.
+ */
+ double Peak::getValueByColName(const std::string & name_in) const
+ {
+ std::string name = name_in;
+ std::transform(name.begin(), name.end(), name.begin(), ::tolower);
+ if (name == "runnumber")
+ return double(this->getRunNumber());
+ else if (name == "detid")
+ return double(this->getDetectorID());
+ else if (name == "h")
+ return this->getH();
+ else if (name == "k")
+ return this->getK();
+ else if (name == "l")
+ return this->getL();
+ else if (name == "wavelength")
+ return this->getWavelength();
+ else if (name == "energy")
+ return this->getInitialEnergy();
+ else if (name == "tof")
+ return this->getTOF();
+ else if (name == "dspacing")
+ return this->getDSpacing();
+ else if (name == "intens")
+ return this->getIntensity();
+ else if (name == "sigint")
+ return this->getSigmaIntensity();
+ else if (name == "bincount")
+ return this->getBinCount();
+ else if (name == "row")
+ return this->getRow();
+ else if (name == "col")
+ return this->getCol();
+ else
+ throw std::runtime_error("Peak::getValueByColName() unknown column or column is not a number: " + name);
+ }
+
+
} // namespace Mantid
} // namespace DataObjects
diff --git a/Code/Mantid/Framework/DataObjects/src/PeakColumn.cpp b/Code/Mantid/Framework/DataObjects/src/PeakColumn.cpp
index 12786c7e670..6f051e3bff8 100644
--- a/Code/Mantid/Framework/DataObjects/src/PeakColumn.cpp
+++ b/Code/Mantid/Framework/DataObjects/src/PeakColumn.cpp
@@ -54,38 +54,14 @@ namespace DataObjects
s << peak.getRunNumber();
else if (m_name == "DetID")
s << peak.getDetectorID();
- else if (m_name == "h")
- s << peak.getH();
- else if (m_name == "k")
- s << peak.getK();
- else if (m_name == "l")
- s << peak.getL();
- else if (m_name == "Wavelength")
- s << peak.getWavelength();
- else if (m_name == "Energy")
- s << peak.getInitialEnergy();
- else if (m_name == "TOF")
- s << peak.getTOF();
- else if (m_name == "DSpacing")
- s << peak.getDSpacing();
- else if (m_name == "Intens")
- s << peak.getIntensity();
- else if (m_name == "SigInt")
- s << peak.getSigmaIntensity();
- else if (m_name == "BinCount")
- s << peak.getBinCount();
else if (m_name == "BankName")
s << peak.getBankName();
- else if (m_name == "Row")
- s << peak.getRow();
- else if (m_name == "Col")
- s << peak.getCol();
else if (m_name == "QLab")
s << peak.getQLabFrame();
else if (m_name == "QSample")
s << peak.getQSampleFrame();
else
- throw std::runtime_error("Unexpected column name");
+ s << peak.getValueByColName(m_name);
}
//-------------------------------------------------------------------------------------
diff --git a/Code/Mantid/Framework/DataObjects/src/PeaksWorkspace.cpp b/Code/Mantid/Framework/DataObjects/src/PeaksWorkspace.cpp
index 10c790d7fb0..57a03c06bac 100644
--- a/Code/Mantid/Framework/DataObjects/src/PeaksWorkspace.cpp
+++ b/Code/Mantid/Framework/DataObjects/src/PeaksWorkspace.cpp
@@ -94,13 +94,85 @@ namespace DataObjects
addColumn( "double", "Intens");
addColumn( "double", "SigInt");
addColumn( "double", "BinCount");
- addColumn( "double", "BankName");
+ addColumn( "str", "BankName");
addColumn( "double", "Row");
addColumn( "double", "Col");
addColumn( "V3D", "QLab");
addColumn( "V3D", "QSample");
}
+ //=====================================================================================
+ //=====================================================================================
+ /** Comparator class for sorting peaks by one or more criteria
+ */
+ class PeakComparator : public std::binary_function<Peak,Peak,bool>
+ {
+ public:
+ std::vector< std::pair<std::string, bool> > & criteria;
+
+ /** Constructor for the comparator for sorting peaks
+ * @param criteria : a vector with a list of pairs: column name, bool;
+ * where bool = true for ascending, false for descending sort.
+ */
+ PeakComparator(std::vector< std::pair<std::string, bool> > & criteria)
+ : criteria(criteria)
+ {
+ }
+
+ /** Compare two peaks using the stored criteria */
+ inline bool operator()(const Peak& a, const Peak& b)
+ {
+ for (size_t i = 0; i < criteria.size(); i++)
+ {
+ std::string & col = criteria[i].first;
+ bool ascending = criteria[i].second;
+ bool lessThan = false;
+ if (col == "BankName")
+ {
+ // If this criterion is equal, move on to the next one
+ std::string valA = a.getBankName();
+ std::string valB = b.getBankName();
+ // Move on to lesser criterion if equal
+ if (valA == valB)
+ continue;
+ lessThan = (valA < valB);
+ }
+ else
+ {
+ // General double comparison
+ double valA = a.getValueByColName(col);
+ double valB = b.getValueByColName(col);
+ // Move on to lesser criterion if equal
+ if (valA == valB)
+ continue;
+ lessThan = (valA < valB);
+ }
+ // Flip the sign of comparison if descending.
+ if (ascending)
+ return lessThan;
+ else
+ return !lessThan;
+
+ }
+ // If you reach here, all criteria were ==; so not <, so return false
+ return false;
+ }
+ };
+
+
+ //---------------------------------------------------------------------------------------------
+ /** Sort the peaks by one or more criteria
+ *
+ * @param criteria : a vector with a list of pairs: column name, bool;
+ * where bool = true for ascending, false for descending sort.
+ * The peaks are sorted by the first criterion first, then the 2nd if equal, etc.
+ */
+ void PeaksWorkspace::sort(std::vector< std::pair<std::string, bool> > & criteria)
+ {
+ PeakComparator comparator(criteria);
+ std::stable_sort(peaks.begin(), peaks.end(), comparator);
+ }
+
//---------------------------------------------------------------------------------------------
/** Destructor */
diff --git a/Code/Mantid/Framework/DataObjects/test/PeakTest.h b/Code/Mantid/Framework/DataObjects/test/PeakTest.h
index 0a5970eda41..4d12b3535be 100644
--- a/Code/Mantid/Framework/DataObjects/test/PeakTest.h
+++ b/Code/Mantid/Framework/DataObjects/test/PeakTest.h
@@ -87,6 +87,21 @@ public:
TS_ASSERT_EQUALS(p.getInstrument(), p2.getInstrument())
}
+ void test_getValueByColName()
+ {
+ Peak p(inst, 10102, 2.0);
+ p.setHKL(1,2,3);
+ p.setRunNumber(1234);
+ TS_ASSERT_EQUALS(p.getValueByColName("Row"), p.getRow());
+ TS_ASSERT_EQUALS(p.getValueByColName("Col"), p.getCol());
+ TS_ASSERT_EQUALS(p.getValueByColName("H"), p.getH());
+ TS_ASSERT_EQUALS(p.getValueByColName("K"), p.getK());
+ TS_ASSERT_EQUALS(p.getValueByColName("L"), p.getL());
+ TS_ASSERT_EQUALS(p.getValueByColName("RunNumber"), p.getRunNumber());
+ TS_ASSERT_EQUALS(p.getValueByColName("DetId"), p.getDetectorID())
+ TS_ASSERT_THROWS_ANYTHING( p.getValueByColName("bankname") );
+ }
+
/** Set the wavelength and see the other "versions" of it get calculated. */
void test_wavelength_conversion()
{
diff --git a/Code/Mantid/Framework/DataObjects/test/PeaksWorkspaceTest.h b/Code/Mantid/Framework/DataObjects/test/PeaksWorkspaceTest.h
index 320570ccc02..ea6532c920c 100644
--- a/Code/Mantid/Framework/DataObjects/test/PeaksWorkspaceTest.h
+++ b/Code/Mantid/Framework/DataObjects/test/PeaksWorkspaceTest.h
@@ -78,6 +78,52 @@ public:
checkPW(pw2.get());
}
+ void test_sort()
+ {
+ PeaksWorkspace_sptr pw(buildPW());
+ Instrument_sptr inst = pw->getInstrument();
+ Peak p0 = pw->getPeak(0); //Peak(inst, 1, 3.0)
+ Peak p1(inst, 1, 4.0);
+ Peak p2(inst, 1, 5.0);
+ Peak p3(inst, 2, 3.0);
+ Peak p4(inst, 3, 3.0);
+ pw->addPeak(p1);
+ pw->addPeak(p2);
+ pw->addPeak(p3);
+ pw->addPeak(p4);
+
+ std::vector< std::pair<std::string, bool> > criteria;
+ // Sort by detector ID then descending wavelength
+ criteria.push_back( std::pair<std::string, bool>("detid", true) );
+ criteria.push_back( std::pair<std::string, bool>("wavelength", false) );
+ pw->sort(criteria);
+ TS_ASSERT_EQUALS( pw->getPeak(0).getDetectorID(), 1);
+ TS_ASSERT_DELTA( pw->getPeak(0).getWavelength(), 5.0, 1e-5);
+ TS_ASSERT_EQUALS( pw->getPeak(1).getDetectorID(), 1);
+ TS_ASSERT_DELTA( pw->getPeak(1).getWavelength(), 4.0, 1e-5);
+ TS_ASSERT_EQUALS( pw->getPeak(2).getDetectorID(), 1);
+ TS_ASSERT_DELTA( pw->getPeak(2).getWavelength(), 3.0, 1e-5);
+ TS_ASSERT_EQUALS( pw->getPeak(3).getDetectorID(), 2);
+ TS_ASSERT_DELTA( pw->getPeak(3).getWavelength(), 3.0, 1e-5);
+
+ // Sort by wavelength ascending then detID
+ criteria.clear();
+ criteria.push_back( std::pair<std::string, bool>("wavelength", true) );
+ criteria.push_back( std::pair<std::string, bool>("detid", true) );
+ pw->sort(criteria);
+ TS_ASSERT_EQUALS( pw->getPeak(0).getDetectorID(), 1);
+ TS_ASSERT_DELTA( pw->getPeak(0).getWavelength(), 3.0, 1e-5);
+ TS_ASSERT_EQUALS( pw->getPeak(1).getDetectorID(), 2);
+ TS_ASSERT_DELTA( pw->getPeak(1).getWavelength(), 3.0, 1e-5);
+ TS_ASSERT_EQUALS( pw->getPeak(2).getDetectorID(), 3);
+ TS_ASSERT_DELTA( pw->getPeak(2).getWavelength(), 3.0, 1e-5);
+ TS_ASSERT_EQUALS( pw->getPeak(3).getDetectorID(), 1);
+ TS_ASSERT_DELTA( pw->getPeak(3).getWavelength(), 4.0, 1e-5);
+ TS_ASSERT_EQUALS( pw->getPeak(4).getDetectorID(), 1);
+ TS_ASSERT_DELTA( pw->getPeak(4).getWavelength(), 5.0, 1e-5);
+
+ }
+
};
diff --git a/Code/Mantid/Framework/MDEvents/src/MDEWFindPeaks.cpp b/Code/Mantid/Framework/MDEvents/src/MDEWFindPeaks.cpp
index 6a56d9685e1..f00eb63c8ea 100644
--- a/Code/Mantid/Framework/MDEvents/src/MDEWFindPeaks.cpp
+++ b/Code/Mantid/Framework/MDEvents/src/MDEWFindPeaks.cpp
@@ -64,22 +64,28 @@ namespace MDEvents
declareProperty(new WorkspaceProperty<IMDEventWorkspace>("InputWorkspace","",Direction::Input),
"An input MDEventWorkspace with at least 3 dimensions.");
- declareProperty(new PropertyWithValue<double>("PeakDistanceThreshold",1.0,Direction::Input),
+ declareProperty(new PropertyWithValue<double>("PeakDistanceThreshold", 0.1, Direction::Input),
"Threshold distance for rejecting peaks that are found to be too close from each other.\n"
- "This should be some multiple of the radius of a peak."
+ "This should be some multiple of the radius of a peak. Default: 0.1."
);
- declareProperty(new PropertyWithValue<int64_t>("MaxPeaks",1000,Direction::Input),
- "Maximum number of peaks to find."
+ declareProperty(new PropertyWithValue<int64_t>("MaxPeaks",500,Direction::Input),
+ "Maximum number of peaks to find. Default: 500."
);
- declareProperty(new PropertyWithValue<double>("DensityThresholdFactor", 1.0, Direction::Input),
- "The overall signal density of the workspace will be multiplied by this factor "
- "to get a threshold signal density below which boxes are NOT considered to be peaks. See the help."
+ declareProperty(new PropertyWithValue<double>("DensityThresholdFactor", 10.0, Direction::Input),
+ "The overall signal density of the workspace will be multiplied by this factor \n"
+ "to get a threshold signal density below which boxes are NOT considered to be peaks. See the help.\n"
+ "Default: 10.0"
);
declareProperty(new WorkspaceProperty<PeaksWorkspace>("OutputWorkspace","",Direction::Output),
"An output PeaksWorkspace with the peaks' found positions.");
+
+ declareProperty("AppendPeaks", false,
+ "If checked, then append the peaks in the output workspace if it exists. \n"
+ "If unchecked, the output workspace is replaced (Default)." );
+
}
@@ -190,7 +196,7 @@ namespace MDEvents
// The box was not rejected for another reason.
if (!badBox)
{
- if (numBoxesFound++ > MaxPeaks)
+ if (numBoxesFound++ >= MaxPeaks)
{
g_log.notice() << "Number of peaks found exceeded the limit of " << MaxPeaks << ". Stopping peak finding." << std::endl;
break;
@@ -250,11 +256,12 @@ namespace MDEvents
*/
void MDEWFindPeaks::exec()
{
- /// Output peaks workspace, create if needed
+ bool AppendPeaks = getProperty("AppendPeaks");
+
+ // Output peaks workspace, create if needed
peakWS = getProperty("OutputWorkspace");
- if (!peakWS)
+ if (!peakWS || !AppendPeaks)
peakWS = PeaksWorkspace_sptr(new PeaksWorkspace());
- setProperty("OutputWorkspace", peakWS);
// The MDEventWorkspace as input
IMDEventWorkspace_sptr inWS = getProperty("InputWorkspace");
@@ -273,6 +280,16 @@ namespace MDEvents
CALL_MDEVENT_FUNCTION3(this->findPeaks, inWS);
delete prog;
+
+ // Do a sort by bank name and then descending bin count (intensity)
+ std::vector< std::pair<std::string, bool> > criteria;
+ criteria.push_back( std::pair<std::string, bool>("BankName", true) );
+ criteria.push_back( std::pair<std::string, bool>("bincount", false) );
+ peakWS->sort(criteria);
+
+ // Save the output
+ setProperty("OutputWorkspace", peakWS);
+
}
diff --git a/Code/Mantid/Framework/MDEvents/test/MDEWFindPeaksTest.h b/Code/Mantid/Framework/MDEvents/test/MDEWFindPeaksTest.h
index e3610e1ba03..81d71f52532 100644
--- a/Code/Mantid/Framework/MDEvents/test/MDEWFindPeaksTest.h
+++ b/Code/Mantid/Framework/MDEvents/test/MDEWFindPeaksTest.h
@@ -69,7 +69,7 @@ public:
TS_ASSERT( alg.isInitialized() )
}
- void test_exec()
+ void do_test(bool deleteWS, int MaxPeaks, int expectedPeaks, bool AppendPeaks = false)
{
// Name of the output workspace.
std::string outWSName("peaksFound");
@@ -89,6 +89,8 @@ public:
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") );
+ TS_ASSERT_THROWS_NOTHING( alg.setProperty("MaxPeaks", int64_t(MaxPeaks)) );
+ TS_ASSERT_THROWS_NOTHING( alg.setProperty("AppendPeaks", AppendPeaks) );
TS_ASSERT_THROWS_NOTHING( alg.execute(); );
TS_ASSERT( alg.isExecuted() );
@@ -99,26 +101,59 @@ public:
if (!ws) return;
// Should find 3 peaks.
- TS_ASSERT_EQUALS( ws->getNumberPeaks(), 3);
- if (ws->getNumberPeaks() != 3) return;
+ TS_ASSERT_EQUALS( ws->getNumberPeaks(), expectedPeaks);
+ if (ws->getNumberPeaks() != expectedPeaks) return;
+ // Stop checking for the AppendPeaks case. This is good enough.
+ if (AppendPeaks) return;
// The order of the peaks found is a little random because it depends on the way the boxes were sorted...
TS_ASSERT_DELTA( ws->getPeak(0).getQLabFrame()[0], -5.0, 0.1);
TS_ASSERT_DELTA( ws->getPeak(0).getQLabFrame()[1], -5.0, 0.1);
TS_ASSERT_DELTA( ws->getPeak(0).getQLabFrame()[2], 5.0, 0.1);
- TS_ASSERT_DELTA( ws->getPeak(1).getQLabFrame()[0], 4.0, 0.1);
- TS_ASSERT_DELTA( ws->getPeak(1).getQLabFrame()[1], 5.0, 0.1);
- TS_ASSERT_DELTA( ws->getPeak(1).getQLabFrame()[2], 6.0, 0.1);
+ if (MaxPeaks > 1)
+ {
+ TS_ASSERT_DELTA( ws->getPeak(1).getQLabFrame()[0], 4.0, 0.1);
+ TS_ASSERT_DELTA( ws->getPeak(1).getQLabFrame()[1], 5.0, 0.1);
+ TS_ASSERT_DELTA( ws->getPeak(1).getQLabFrame()[2], 6.0, 0.1);
+
+ TS_ASSERT_DELTA( ws->getPeak(2).getQLabFrame()[0], 1.0, 0.1);
+ TS_ASSERT_DELTA( ws->getPeak(2).getQLabFrame()[1], 2.0, 0.1);
+ TS_ASSERT_DELTA( ws->getPeak(2).getQLabFrame()[2], 3.0, 0.1);
+ }
+
+ if (deleteWS)
+ {
+ // Remove workspace from the data service.
+ AnalysisDataService::Instance().remove(outWSName);
+ }
+ AnalysisDataService::Instance().remove("MDEWS");
+ }
- TS_ASSERT_DELTA( ws->getPeak(2).getQLabFrame()[0], 1.0, 0.1);
- TS_ASSERT_DELTA( ws->getPeak(2).getQLabFrame()[1], 2.0, 0.1);
- TS_ASSERT_DELTA( ws->getPeak(2).getQLabFrame()[2], 3.0, 0.1);
-
- // Remove workspace from the data service.
- AnalysisDataService::Instance().remove(outWSName);
+ /** Running the algo twice with same output workspace = replace the output, don't append */
+ void test_exec_twice_replaces_workspace()
+ {
+ do_test(false, 100, 3);
+ do_test(true, 100, 3);
}
-
+
+ /** Run twice and append to the peaks workspace*/
+ void test_exec_AppendPeaks()
+ {
+ do_test(false, 100, 3);
+ do_test(true, 100, 6, true /* Append */ );
+ }
+
+ void test_exec()
+ {
+ do_test(true, 100, 3);
+ }
+
+ void test_exec_withMaxPeaks()
+ {
+ do_test(true, 1, 1);
+ }
+
};
--
GitLab