Commit 07fed966 authored by Janik Zikovsky's avatar Janik Zikovsky
Browse files

Refs #3590: MDEWFindPeaks: AppendPeaks option, replace output if not checked....

Refs #3590: MDEWFindPeaks: AppendPeaks option, replace output if not checked. Also sorting output by bankname and intensity.
parent 375ace2f
......@@ -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;
......
......@@ -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
*/
......
......@@ -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
......@@ -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);
}
//-------------------------------------------------------------------------------------
......
......@@ -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 */
......
......@@ -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()
{
......
......@@ -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);
}
};
......
......@@ -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);
}
......
......@@ -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);
}
};
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment