Commit 84587983 authored by Lynch, Vickie's avatar Lynch, Vickie
Browse files

Refs #8241 new versions

parent 8b1001ad
......@@ -113,13 +113,6 @@ namespace Crystal
declareProperty(new WorkspaceProperty<PeaksWorkspace>("PeaksWorkspace","",Direction::Input),
"The list of peaks to integrate, matching the InputWorkspace.");
std::vector<std::string> propOptions;
propOptions.push_back("Q (lab frame)");
propOptions.push_back("Q (sample frame)");
propOptions.push_back("HKL");
declareProperty("CoordinatesToUse", "Q (lab frame)", boost::make_shared<StringListValidator>(propOptions),
"Deprecated: algorithm uses the InputWorkspace's coordinates.");
declareProperty("RadiusStart", 0.0, "Radius at which to start integrating." );
declareProperty("RadiusEnd", 1.0, "Radius at which to stop integrating." );
declareProperty("NumSteps", 10, "Number of steps, between start and end, to calculate radius." );
......@@ -235,7 +228,6 @@ namespace Crystal
IAlgorithm_sptr alg = this->createChildAlgorithm("IntegratePeaksMD", progStep*double(step), progStep*double(step+1), false);
alg->setProperty("InputWorkspace", inWS);
alg->setProperty("PeaksWorkspace", peaksWS);
alg->setPropertyValue("CoordinatesToUse", this->getPropertyValue("CoordinatesToUse"));
alg->setProperty("PeakRadius", radius);
alg->setProperty("BackgroundOuterRadius", OuterRadius);
alg->setProperty("BackgroundInnerRadius", InnerRadius);
......
......@@ -82,7 +82,6 @@ public:
TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("InputWorkspace", "PeakIntensityVsRadiusTest_MDEWS") );
TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("PeaksWorkspace", "PeakIntensityVsRadiusTest_peaks") );
TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("OutputWorkspace", "PeakIntensityVsRadiusTest_OutputWS") );
TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("CoordinatesToUse", "HKL") );
TS_ASSERT_THROWS_NOTHING( alg.setProperty("RadiusStart", 0.0) );
TS_ASSERT_THROWS_NOTHING( alg.setProperty("RadiusEnd", 1.5) );
TS_ASSERT_THROWS_NOTHING( alg.setProperty("NumSteps", 16) );
......@@ -121,7 +120,6 @@ public:
TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("InputWorkspace", "PeakIntensityVsRadiusTest_MDEWS") );
TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("PeaksWorkspace", "PeakIntensityVsRadiusTest_peaks") );
TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("OutputWorkspace", outWSName) );
TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("CoordinatesToUse", "HKL") );
TS_ASSERT_THROWS_NOTHING( alg.setProperty("RadiusStart", 0.0) );
TS_ASSERT_THROWS_NOTHING( alg.setProperty("RadiusEnd", 1.5) );
TS_ASSERT_THROWS_NOTHING( alg.setProperty("NumSteps", 16) );
......
......@@ -8,6 +8,7 @@ set ( SRC_FILES
src/BinaryOperationMD.cpp
src/BooleanBinaryOperationMD.cpp
src/CentroidPeaksMD.cpp
src/CentroidPeaksMD2.cpp
src/CloneMDWorkspace.cpp
src/CompareMDWorkspaces.cpp
src/ConvertToDetectorFaceMD.cpp
......@@ -25,6 +26,7 @@ set ( SRC_FILES
src/GreaterThanMD.cpp
src/IDynamicRebinning.cpp
src/IntegratePeaksMD.cpp
src/IntegratePeaksMD2.cpp
src/InvalidParameter.cpp
src/InvalidParameterParser.cpp
src/LessThanMD.cpp
......@@ -78,6 +80,7 @@ set ( INC_FILES
inc/MantidMDAlgorithms/BinaryOperationMD.h
inc/MantidMDAlgorithms/BooleanBinaryOperationMD.h
inc/MantidMDAlgorithms/CentroidPeaksMD.h
inc/MantidMDAlgorithms/CentroidPeaksMD2.h
inc/MantidMDAlgorithms/CloneMDWorkspace.h
inc/MantidMDAlgorithms/CompareMDWorkspaces.h
inc/MantidMDAlgorithms/ConvertToDetectorFaceMD.h
......@@ -97,6 +100,7 @@ set ( INC_FILES
inc/MantidMDAlgorithms/GSLFunctions.h
inc/MantidMDAlgorithms/IDynamicRebinning.h
inc/MantidMDAlgorithms/IntegratePeaksMD.h
inc/MantidMDAlgorithms/IntegratePeaksMD2.h
inc/MantidMDAlgorithms/InvalidParameter.h
inc/MantidMDAlgorithms/InvalidParameterParser.h
inc/MantidMDAlgorithms/LessThanMD.h
......@@ -150,6 +154,7 @@ set ( TEST_FILES
BooleanBinaryOperationMDTest.h
CachedExperimentInfoTest.h
CentroidPeaksMDTest.h
CentroidPeaksMD2Test.h
CloneMDWorkspaceTest.h
CompareMDWorkspacesTest.h
ConvertEventsToMDTest.h
......@@ -171,6 +176,7 @@ set ( TEST_FILES
ForegroundModelTest.h
GreaterThanMDTest.h
IntegratePeaksMDTest.h
IntegratePeaksMD2Test.h
InvalidParameterParserTest.h
InvalidParameterTest.h
LessThanMDTest.h
......
#ifndef MANTID_MDALGORITHMS_CENTROIDPEAKSMD2_H_
#define MANTID_MDALGORITHMS_CENTROIDPEAKSMD2_H_
#include "MantidKernel/System.h"
#include "MantidAPI/Algorithm.h"
#include "MantidAPI/IMDEventWorkspace.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidMDEvents/MDEventWorkspace.h"
namespace Mantid
{
namespace MDAlgorithms
{
/** Find the centroid of single-crystal peaks in a MDEventWorkspace, in order to refine their positions.
*
* @author Janik Zikovsky
* @date 2011-06-01
*/
class DLLExport CentroidPeaksMD2 : public API::Algorithm
{
public:
CentroidPeaksMD2();
~CentroidPeaksMD2();
/// Algorithm's name for identification
virtual const std::string name() const { return "CentroidPeaksMD";};
/// Algorithm's version for identification
virtual int version() const { return 2;};
/// Algorithm's category for identification
virtual const std::string category() const { return "MDAlgorithms";}
private:
/// Sets documentation strings for this algorithm
virtual void initDocs();
/// Initialise the properties
void init();
/// Run the algorithm
void exec();
template<typename MDE, size_t nd>
void integrate(typename MDEvents::MDEventWorkspace<MDE, nd>::sptr ws);
/// Input MDEventWorkspace
Mantid::API::IMDEventWorkspace_sptr inWS;
};
} // namespace Mantid
} // namespace MDEvents
#endif /* MANTID_MDEVENTS_CENTROIDPEAKSMD2_H_ */
......@@ -7,6 +7,7 @@ namespace Mantid
namespace MDAlgorithms
{
double f_eval (double x, void *params);
double f_eval2 (double x, void *params);
} // namespace MDAlgorithms
} // namespace Mantid
......
#ifndef MANTID_MDALGORITHMS_INTEGRATEPEAKSMD_H_
#define MANTID_MDALGORITHMS_INTEGRATEPEAKSMD_H_
#include "MantidAPI/Algorithm.h"
#include "MantidAPI/IMDEventWorkspace.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidKernel/System.h"
#include "MantidMDEvents/MDEventWorkspace.h"
#include "MantidAPI/CompositeFunction.h"
namespace Mantid
{
namespace MDAlgorithms
{
/** Integrate single-crystal peaks in reciprocal-space.
*
* @author Janik Zikovsky
* @date 2011-04-13 18:11:53.496539
*/
class DLLExport IntegratePeaksMD2 : public API::Algorithm
{
public:
IntegratePeaksMD2();
~IntegratePeaksMD2();
/// Algorithm's name for identification
virtual const std::string name() const { return "IntegratePeaksMD";};
/// Algorithm's version for identification
virtual int version() const { return 2;};
/// Algorithm's category for identification
virtual const std::string category() const { return "MDAlgorithms";}
private:
/// Sets documentation strings for this algorithm
virtual void initDocs();
/// Initialise the properties
void init();
/// Run the algorithm
void exec();
template<typename MDE, size_t nd>
void integrate(typename MDEvents::MDEventWorkspace<MDE, nd>::sptr ws);
/// Input MDEventWorkspace
Mantid::API::IMDEventWorkspace_sptr inWS;
/// Calculate if this Q is on a detector
bool detectorQ(Mantid::Kernel::V3D QLabFrame, double PeakRadius);
/// Instrument reference
Geometry::Instrument_const_sptr inst;
};
} // namespace Mantid
} // namespace MDEvents
#endif /* MANTID_MDALGORITHMS_INTEGRATEPEAKSMD_H_ */
......@@ -67,7 +67,7 @@ namespace MDAlgorithms
propOptions.push_back("Q (sample frame)");
propOptions.push_back("HKL");
declareProperty("CoordinatesToUse", "HKL",boost::make_shared<StringListValidator>(propOptions),
"Deprecated: algorithm uses the InputWorkspace's coordinates.");
"Ignored: algorithm uses the InputWorkspace's coordinates.");
declareProperty(new PropertyWithValue<double>("PeakRadius",1.0,Direction::Input),
"Fixed radius around each peak position in which to calculate the centroid.");
......
/*WIKI*
This algorithm starts with a PeaksWorkspace containing the expected positions of peaks in reciprocal space. It calculates the centroid of the peak by calculating the average of the coordinates of all events within a given radius of the peak, weighted by the weight (signal) of the event.
*WIKI*/
#include "MantidKernel/System.h"
#include "MantidKernel/ListValidator.h"
#include "MantidAPI/IMDEventWorkspace.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidMDEvents/CoordTransformDistance.h"
#include "MantidMDEvents/MDEventFactory.h"
#include "MantidMDAlgorithms/IntegratePeaksMD.h"
#include "MantidMDAlgorithms/CentroidPeaksMD2.h"
using Mantid::DataObjects::PeaksWorkspace;
namespace Mantid
{
namespace MDAlgorithms
{
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(CentroidPeaksMD2)
using namespace Mantid::API;
using namespace Mantid::DataObjects;
using namespace Mantid::Geometry;
using namespace Mantid::Kernel;
using namespace Mantid::MDEvents;
//----------------------------------------------------------------------------------------------
/** Constructor
*/
CentroidPeaksMD2::CentroidPeaksMD2()
{
}
//----------------------------------------------------------------------------------------------
/** Destructor
*/
CentroidPeaksMD2::~CentroidPeaksMD2()
{
}
//----------------------------------------------------------------------------------------------
/// Sets documentation strings for this algorithm
void CentroidPeaksMD2::initDocs()
{
this->setWikiSummary("Find the centroid of single-crystal peaks in a MDEventWorkspace, in order to refine their positions.");
this->setOptionalMessage("Find the centroid of single-crystal peaks in a MDEventWorkspace, in order to refine their positions.");
}
//----------------------------------------------------------------------------------------------
/** Initialize the algorithm's properties.
*/
void CentroidPeaksMD2::init()
{
declareProperty(new WorkspaceProperty<IMDEventWorkspace>("InputWorkspace","",Direction::Input), "An input MDEventWorkspace.");
std::vector<std::string> propOptions;
propOptions.push_back("Q (lab frame)");
propOptions.push_back("Q (sample frame)");
propOptions.push_back("HKL");
declareProperty(new PropertyWithValue<double>("PeakRadius",1.0,Direction::Input),
"Fixed radius around each peak position in which to calculate the centroid.");
declareProperty(new WorkspaceProperty<PeaksWorkspace>("PeaksWorkspace","",Direction::Input),
"A PeaksWorkspace containing the peaks to centroid.");
declareProperty(new WorkspaceProperty<PeaksWorkspace>("OutputWorkspace","",Direction::Output),
"The output PeaksWorkspace will be a copy of the input PeaksWorkspace "
"with the peaks' positions modified by the new found centroids.");
}
//----------------------------------------------------------------------------------------------
/** Integrate the peaks of the workspace using parameters saved in the algorithm class
* @param ws :: MDEventWorkspace to integrate
*/
template<typename MDE, size_t nd>
void CentroidPeaksMD2::integrate(typename MDEventWorkspace<MDE, nd>::sptr ws)
{
if (nd != 3)
throw std::invalid_argument("For now, we expect the input MDEventWorkspace to have 3 dimensions only.");
/// Peak workspace to centroid
Mantid::DataObjects::PeaksWorkspace_sptr inPeakWS = getProperty("PeaksWorkspace");
/// Output peaks workspace, create if needed
Mantid::DataObjects::PeaksWorkspace_sptr peakWS = getProperty("OutputWorkspace");
if (peakWS != inPeakWS)
peakWS = inPeakWS->clone();
int CoordinatesToUse = ws->getSpecialCoordinateSystem();
/// Radius to use around peaks
double PeakRadius = getProperty("PeakRadius");
// cppcheck-suppress syntaxError
PRAGMA_OMP(parallel for schedule(dynamic, 10) )
for (int i=0; i < int(peakWS->getNumberPeaks()); ++i)
{
// Get a direct ref to that peak.
IPeak & p = peakWS->getPeak(i);
double detectorDistance = p.getL2();
// Get the peak center as a position in the dimensions of the workspace
V3D pos;
if (CoordinatesToUse == 1) //"Q (lab frame)"
pos = p.getQLabFrame();
else if (CoordinatesToUse == 2) //"Q (sample frame)"
pos = p.getQSampleFrame();
else if (CoordinatesToUse == 3) //"HKL"
pos = p.getHKL();
// Build the sphere transformation
bool dimensionsUsed[nd];
coord_t center[nd];
for (size_t d=0; d<nd; ++d)
{
dimensionsUsed[d] = true; // Use all dimensions
center[d] = static_cast<coord_t>(pos[d]);
}
CoordTransformDistance sphere(nd, center, dimensionsUsed);
// Initialize the centroid to 0.0
signal_t signal = 0;
coord_t centroid[nd];
for (size_t d=0; d<nd; d++)
centroid[d] = 0.0;
// Perform centroid
ws->getBox()->centroidSphere(sphere, static_cast<coord_t>(PeakRadius*PeakRadius), centroid, signal);
// Normalize by signal
if (signal != 0.0)
{
for (size_t d=0; d<nd; d++)
centroid[d] /= static_cast<coord_t>(signal);
V3D vecCentroid(centroid[0], centroid[1], centroid[2]);
// Save it back in the peak object, in the dimension specified.
if (CoordinatesToUse == 1) //"Q (lab frame)"
{
p.setQLabFrame( vecCentroid, detectorDistance);
p.findDetector();
}
else if (CoordinatesToUse == 2) //"Q (sample frame)"
{
p.setQSampleFrame( vecCentroid, detectorDistance);
p.findDetector();
}
else if (CoordinatesToUse == 3) //"HKL"
{
p.setHKL( vecCentroid );
}
g_log.information() << "Peak " << i << " at " << pos << ": signal "
<< signal << ", centroid " << vecCentroid
<< " in " << CoordinatesToUse
<< std::endl;
}
else
{
g_log.information() << "Peak " << i << " at " << pos << " had no signal, and could not be centroided."
<< std::endl;
}
}
// Save the output
setProperty("OutputWorkspace", peakWS);
}
//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
*/
void CentroidPeaksMD2::exec()
{
inWS = getProperty("InputWorkspace");
CALL_MDEVENT_FUNCTION3(this->integrate, inWS);
}
} // namespace Mantid
} // namespace MDEvents
......@@ -139,7 +139,7 @@ namespace MDAlgorithms
propOptions.push_back("Q (sample frame)");
propOptions.push_back("HKL");
declareProperty("CoordinatesToUse", "Q (lab frame)",boost::make_shared<StringListValidator>(propOptions),
"Deprecated: algorithm uses the InputWorkspace's coordinates.");
"Ignored: algorithm uses the InputWorkspace's coordinates.");
declareProperty(new PropertyWithValue<double>("PeakRadius",1.0,Direction::Input),
"Fixed radius around each peak position in which to integrate (in the same units as the workspace).");
......
This diff is collapsed.
#ifndef MANTID_MDEVENTS_MDCENTROIDPEAKS2TEST_H_
#define MANTID_MDEVENTS_MDCENTROIDPEAKS2TEST_H_
#include "MantidMDAlgorithms/CentroidPeaksMD2.h"
#include "MantidAPI/AnalysisDataService.h"
#include "MantidAPI/IMDEventWorkspace.h"
#include "MantidAPI/FrameworkManager.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidGeometry/MDGeometry/MDHistoDimension.h"
#include "MantidKernel/System.h"
#include "MantidKernel/Timer.h"
#include "MantidMDEvents/MDEventFactory.h"
#include "MantidTestHelpers/ComponentCreationHelper.h"
#include <boost/math/distributions/normal.hpp>
#include <boost/math/special_functions/fpclassify.hpp>
#include <boost/math/special_functions/pow.hpp>
#include <boost/random/linear_congruential.hpp>
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/uniform_int.hpp>
#include <boost/random/uniform_real.hpp>
#include <boost/random/variate_generator.hpp>
#include <cxxtest/TestSuite.h>
#include <iomanip>
#include <iostream>
using Mantid::API::AnalysisDataService;
using Mantid::Geometry::MDHistoDimension;
using namespace Mantid::API;
using namespace Mantid::DataObjects;
using namespace Mantid::Geometry;
using namespace Mantid::MDEvents;
using namespace Mantid::MDAlgorithms;
using Mantid::Kernel::V3D;
class CentroidPeaksMD2Test : public CxxTest::TestSuite
{
public:
void test_Init()
{
CentroidPeaksMD2 alg;
TS_ASSERT_THROWS_NOTHING( alg.initialize() )
TS_ASSERT( alg.isInitialized() )
}
//-------------------------------------------------------------------------------
/** Create the (blank) MDEW */
static void createMDEW()
{
// ---- Start with empty MDEW ----
FrameworkManager::Instance().exec("CreateMDWorkspace", 14,
"Dimensions", "3",
"Extents", "-10,10,-10,10,-10,10",
"Names", "h,k,l",
"Units", "-,-,-",
"SplitInto", "5",
"MaxRecursionDepth", "2",
"OutputWorkspace", "CentroidPeaksMD2Test_MDEWS");
}
//-------------------------------------------------------------------------------
/** Add a fake "peak"*/
static void addPeak(size_t num, double x, double y, double z, double radius)
{
std::ostringstream mess;
mess << num << ", " << x << ", " << y << ", " << z << ", " << radius;
FrameworkManager::Instance().exec("FakeMDEventData", 6,
"InputWorkspace", "CentroidPeaksMD2Test_MDEWS",
"PeakParams", mess.str().c_str(),
"RandomSeed", "1234");
}
//-------------------------------------------------------------------------------
/** Run the CentroidPeaksMD2 with the given peak radius param */
void doRun( V3D startPos, double PeakRadius, V3D expectedResult, std::string message,
std::string OutputWorkspace = "CentroidPeaksMD2Test_Peaks")
{
// Make a fake instrument - doesn't matter, we won't use it really
Instrument_sptr inst = ComponentCreationHelper::createTestInstrumentCylindrical(5);
// --- Make a fake PeaksWorkspace in the given coordinate space ---
PeaksWorkspace_sptr peakWS(new PeaksWorkspace());
Peak pIn(inst, 1, 1.0, startPos );
if (CoordinatesToUse == "Q (lab frame)")
pIn.setQLabFrame(startPos);
else if (CoordinatesToUse == "Q (sample frame)")
pIn.setQSampleFrame(startPos);
else if (CoordinatesToUse == "HKL")
pIn.setHKL(startPos);
peakWS->addPeak( pIn );
TS_ASSERT_EQUALS( peakWS->getPeak(0).getIntensity(), 0.0);
AnalysisDataService::Instance().addOrReplace("CentroidPeaksMD2Test_Peaks", peakWS);
CentroidPeaksMD2 alg;
TS_ASSERT_THROWS_NOTHING( alg.initialize() )
TS_ASSERT( alg.isInitialized() )
TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("InputWorkspace", "CentroidPeaksMD2Test_MDEWS" ) );
TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("PeaksWorkspace", "CentroidPeaksMD2Test_Peaks" ) );
TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("OutputWorkspace", OutputWorkspace ) );
TS_ASSERT_THROWS_NOTHING( alg.setProperty("PeakRadius", PeakRadius ) );
TS_ASSERT_THROWS_NOTHING( alg.execute() );
TS_ASSERT( alg.isExecuted() );
peakWS = boost::dynamic_pointer_cast<PeaksWorkspace>(
AnalysisDataService::Instance().retrieve(OutputWorkspace));
TS_ASSERT(peakWS);
if (!peakWS) return;
// Compare the result to the expectation
V3D result;
IPeak & p = peakWS->getPeak(0);
if (CoordinatesToUse == "Q (lab frame)")
result = p.getQLabFrame();
else if (CoordinatesToUse == "Q (sample frame)")
{
std::cerr << p.getGoniometerMatrix() << std::endl;
result = p.getQSampleFrame();
}
else if (CoordinatesToUse == "HKL")
result = p.getHKL();
for (size_t i=0; i<3; i++)
TSM_ASSERT_DELTA( message, result[i], expectedResult[i], 0.05);
AnalysisDataService::Instance().remove("CentroidPeaksMD2Test_Peaks");
}
//-------------------------------------------------------------------------------
/** Full test using faked-out peak data */
void do_test_exec()
{
// --- Fake workspace with 3 peaks ------
createMDEW();
addPeak(1000, 0,0.,0., 1.0);
addPeak(1000, 2.,3.,4., 0.5);
addPeak(1000, 6.,6.,6., 2.0);
MDEventWorkspace3Lean::sptr mdews = AnalysisDataService::Instance().retrieveWS<MDEventWorkspace3Lean>("CentroidPeaksMD2Test_MDEWS");
TS_ASSERT_EQUALS( mdews->getNPoints(), 3000);
TS_ASSERT_DELTA( mdews->getBox()->getSignal(), 3000.0, 1e-2);
if (CoordinatesToUse == "HKL")
{
mdews->setCoordinateSystem(Mantid::API::HKL);
doRun(V3D( 0.,0.,0.), 1.0, V3D( 0.,0.,0.), "Start at the center, get the center");
doRun(V3D( 0.2,0.2,0.2), 1.8, V3D( 0.,0.,0.), "Somewhat off center");
}
else if (CoordinatesToUse == "Q (lab frame)")
{
mdews->setCoordinateSystem(Mantid::API::QLab);
}
else if (CoordinatesToUse == "Q (sample frame)")