diff --git a/Code/Mantid/Framework/SINQ/CMakeLists.txt b/Code/Mantid/Framework/SINQ/CMakeLists.txt
index 9a4e5197d75a..cec1ba898692 100644
--- a/Code/Mantid/Framework/SINQ/CMakeLists.txt
+++ b/Code/Mantid/Framework/SINQ/CMakeLists.txt
@@ -13,6 +13,7 @@ set ( SRC_FILES
src/PoldiPeakSearch.cpp
src/PoldiRemoveDeadWires.cpp
src/PoldiUtilities/MillerIndices.cpp
+ src/PoldiUtilities/PeakFunctionIntegrator.cpp
src/PoldiUtilities/PoldiAutoCorrelationCore.cpp
src/PoldiUtilities/PoldiBasicChopper.cpp
src/PoldiUtilities/PoldiChopperFactory.cpp
@@ -52,6 +53,7 @@ set ( INC_FILES
inc/MantidSINQ/PoldiRemoveDeadWires.h
inc/MantidSINQ/PoldiUtilities/MillerIndices.h
inc/MantidSINQ/PoldiUtilities/MillerIndicesIO.h
+ inc/MantidSINQ/PoldiUtilities/PeakFunctionIntegrator.h
inc/MantidSINQ/PoldiUtilities/PoldiAbstractChopper.h
inc/MantidSINQ/PoldiUtilities/PoldiAbstractDetector.h
inc/MantidSINQ/PoldiUtilities/PoldiAutoCorrelationCore.h
@@ -83,6 +85,7 @@ set ( TEST_FILES
MDHistoToWorkspace2DTest.h
MillerIndicesIOTest.h
MillerIndicesTest.h
+ PeakFunctionIntegratorTest.h
PoldiAutoCorrelationCoreTest.h
PoldiBasicChopperTest.h
PoldiChopperFactoryTest.h
diff --git a/Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiUtilities/PeakFunctionIntegrator.h b/Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiUtilities/PeakFunctionIntegrator.h
new file mode 100644
index 000000000000..c9c868ee91c4
--- /dev/null
+++ b/Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiUtilities/PeakFunctionIntegrator.h
@@ -0,0 +1,79 @@
+#ifndef PEAKFUNCTIONINTEGRATOR_H
+#define PEAKFUNCTIONINTEGRATOR_H
+
+#include "MantidSINQ/DllConfig.h"
+#include "MantidAPI/IPeakFunction.h"
+#include "gsl/gsl_integration.h"
+
+namespace Mantid {
+namespace Poldi {
+
+/** PeakFunctionIntegrator :
+ *
+ General integration of peaks (in the form of IPeakFunction) by wrapping the
+ corresponding GSL-functions. Integration with infinity limits is supported.
+
+ @author Michael Wedel, Paul Scherrer Institut - SINQ
+ @date 24/04/2014
+
+ Copyright © 2014 PSI-MSS
+
+ This file is part of Mantid.
+
+ Mantid is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ Mantid is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see .
+
+ File change history is stored at:
+ Code Documentation is available at:
+*/
+
+using namespace Mantid::API;
+
+struct MANTID_SINQ_DLL IntegrationResult {
+ double result;
+ double error;
+ size_t intervals;
+
+ int errorCode;
+ bool success;
+};
+
+class MANTID_SINQ_DLL PeakFunctionIntegrator {
+public:
+ PeakFunctionIntegrator(double requiredRelativePrecision = 1e-8);
+ virtual ~PeakFunctionIntegrator();
+
+ void setRequiredRelativePrecision(double newPrecision);
+ double requiredRelativePrecision() const;
+
+ IntegrationResult integrateInfinity(IPeakFunction_const_sptr peakFunction) const;
+ IntegrationResult integratePositiveInfinity(IPeakFunction_const_sptr peakFunction, double lowerLimit) const;
+ IntegrationResult integrateNegativeInfinity(IPeakFunction_const_sptr peakFunction, double upperLimit) const;
+
+ IntegrationResult integrate(IPeakFunction_const_sptr peakFunction, double lowerLimit, double upperLimit) const;
+
+protected:
+ gsl_function getGSLFunction(IPeakFunction_const_sptr peakFunction) const;
+ void throwIfInvalid(IPeakFunction_const_sptr peakFunction) const;
+
+ gsl_integration_workspace *m_integrationWorkspace;
+
+ double m_relativePrecision;
+};
+
+double gsl_peak_wrapper(double x, void *parameters);
+
+}
+}
+
+#endif // PEAKFUNCTIONINTEGRATOR_H
diff --git a/Code/Mantid/Framework/SINQ/src/PoldiUtilities/PeakFunctionIntegrator.cpp b/Code/Mantid/Framework/SINQ/src/PoldiUtilities/PeakFunctionIntegrator.cpp
new file mode 100644
index 000000000000..d557a5a4a93e
--- /dev/null
+++ b/Code/Mantid/Framework/SINQ/src/PoldiUtilities/PeakFunctionIntegrator.cpp
@@ -0,0 +1,173 @@
+#include "MantidSINQ/PoldiUtilities/PeakFunctionIntegrator.h"
+
+#include "MantidAPI/FunctionDomain1D.h"
+#include "gsl/gsl_errno.h"
+#include
+#include
+
+namespace Mantid {
+namespace Poldi {
+
+/** Constructor with required relative precision argument. The default is 1e-8.
+ * See also PeakFunctionIntegrator::setRequiredRelativePrecision.
+ *
+ * @param requiredRelativePrecision :: Desired relative precision of the integral estimations.
+ */
+PeakFunctionIntegrator::PeakFunctionIntegrator(double requiredRelativePrecision) :
+ m_integrationWorkspace(gsl_integration_workspace_alloc(1000)),
+ m_relativePrecision(requiredRelativePrecision)
+{
+ /* Error handling is disabled, so error-codes are associated to the integration result
+ * and have to be checked.
+ */
+ gsl_set_error_handler_off();
+}
+
+PeakFunctionIntegrator::~PeakFunctionIntegrator()
+{
+ gsl_integration_workspace_free(m_integrationWorkspace);
+}
+
+/** This method sets the desired numerical relative precision that's passed on to the
+ * GSL integration-routines.
+ *
+ * @param newPrecision :: Desired relative precision for integrations.
+ */
+void PeakFunctionIntegrator::setRequiredRelativePrecision(double newPrecision)
+{
+ m_relativePrecision = newPrecision;
+}
+
+/** Returns the currently set precision
+ */
+double PeakFunctionIntegrator::requiredRelativePrecision() const
+{
+ return m_relativePrecision;
+}
+
+/** Integration of peak function on the interval [-Inf, +Inf]. Internally, gsl_integration_qagi is used
+ * for this. If a default constructed IPeakFunction_const_sptr is passed to the function, std::invalid_argument is thrown.
+ * The results are returned as IntegrationResult-struct, which contains the approximation of the integral along
+ * with other information such as an error estimate (absolute).
+ *
+ * @param peakFunction :: Peak function to integrate.
+ */
+IntegrationResult PeakFunctionIntegrator::integrateInfinity(IPeakFunction_const_sptr peakFunction) const
+{
+ throwIfInvalid(peakFunction);
+
+ IntegrationResult result;
+
+ gsl_function f = getGSLFunction(peakFunction);
+
+ result.errorCode = gsl_integration_qagi(&f, 0, m_relativePrecision, 1000, m_integrationWorkspace, &result.result, &result.error);
+ result.success = (result.errorCode == GSL_SUCCESS);
+ result.intervals = m_integrationWorkspace->size;
+
+ return result;
+}
+
+/** Integration of peak function on the interval [a, +Inf]. Internally, gsl_integration_qagiu is used
+ * for this. If a default constructed IPeakFunction_const_sptr is passed to the function, std::invalid_argument is thrown.
+ *
+ * @param peakFunction :: Peak function to integrate.
+ */
+IntegrationResult PeakFunctionIntegrator::integratePositiveInfinity(IPeakFunction_const_sptr peakFunction, double lowerLimit) const
+{
+ throwIfInvalid(peakFunction);
+
+ IntegrationResult result;
+
+ gsl_function f = getGSLFunction(peakFunction);
+
+ result.errorCode = gsl_integration_qagiu(&f, lowerLimit, 0, m_relativePrecision, 1000, m_integrationWorkspace, &result.result, &result.error);
+ result.success = (result.errorCode == GSL_SUCCESS);
+ result.intervals = m_integrationWorkspace->size;
+
+ return result;
+}
+
+/** Integration of peak function on the interval [-Inf, b]. Internally, gsl_integration_qagil is used
+ * for this. If a default constructed IPeakFunction_const_sptr is passed to the function, std::invalid_argument is thrown.
+ *
+ * @param peakFunction :: Peak function to integrate.
+ */
+IntegrationResult PeakFunctionIntegrator::integrateNegativeInfinity(IPeakFunction_const_sptr peakFunction, double upperLimit) const
+{
+ throwIfInvalid(peakFunction);
+
+ IntegrationResult result;
+
+ gsl_function f = getGSLFunction(peakFunction);
+
+ result.errorCode = gsl_integration_qagil(&f, upperLimit, 0, m_relativePrecision, 1000, m_integrationWorkspace, &result.result, &result.error);
+ result.success = (result.errorCode == GSL_SUCCESS);
+ result.intervals = m_integrationWorkspace->size;
+
+ return result;
+}
+
+/** Integration of peak function on the interval [a, b]. Internally, gsl_integration_qags is used
+ * for this. If a default constructed IPeakFunction_const_sptr is passed to the function, std::invalid_argument is thrown.
+ *
+ * @param peakFunction :: Peak function to integrate.
+ */
+IntegrationResult PeakFunctionIntegrator::integrate(IPeakFunction_const_sptr peakFunction, double lowerLimit, double upperLimit) const
+{
+ throwIfInvalid(peakFunction);
+
+ IntegrationResult result;
+
+ gsl_function f = getGSLFunction(peakFunction);
+
+ result.errorCode = gsl_integration_qags(&f, lowerLimit, upperLimit, 0, m_relativePrecision, 1000, m_integrationWorkspace, &result.result, &result.error);
+ result.success = (result.errorCode == GSL_SUCCESS);
+ result.intervals = m_integrationWorkspace->size;
+
+ return result;
+}
+
+/** Method that wraps an IPeakFunction for use with GSL functions.
+ *
+ * @param peakFunction :: Peak function to wrap.
+ */
+gsl_function PeakFunctionIntegrator::getGSLFunction(IPeakFunction_const_sptr peakFunction) const
+{
+ gsl_function f;
+ f.function = &Mantid::Poldi::gsl_peak_wrapper;
+ f.params = &peakFunction;
+
+ return f;
+}
+
+void PeakFunctionIntegrator::throwIfInvalid(IPeakFunction_const_sptr peakFunction) const
+{
+ if(!peakFunction) {
+ throw std::invalid_argument("Can not integrate NULL-function.");
+ }
+}
+
+double gsl_peak_wrapper(double x, void *parameters)
+{
+ IPeakFunction_const_sptr peakFunction = *(IPeakFunction_const_sptr *) parameters;
+
+ double y;
+
+ /* For the integration to work properly, functionLocal has to be used instead
+ * of the more general function-method due to the fact that the overriden function-method
+ * in IPeakFunction cuts off at some point. For slowly decaying peak functions
+ * such as Lorentzians, this is introduces large deviations for integrations
+ * from -Inf to +Inf.
+ */
+ peakFunction->functionLocal(&y, &x, 1);
+
+ return y;
+}
+
+
+
+}
+
+
+
+}
diff --git a/Code/Mantid/Framework/SINQ/test/PeakFunctionIntegratorTest.h b/Code/Mantid/Framework/SINQ/test/PeakFunctionIntegratorTest.h
new file mode 100644
index 000000000000..c01960d54b29
--- /dev/null
+++ b/Code/Mantid/Framework/SINQ/test/PeakFunctionIntegratorTest.h
@@ -0,0 +1,187 @@
+#ifndef PEAKFUNCTIONINTEGRATORTEST_H
+#define PEAKFUNCTIONINTEGRATORTEST_H
+
+#include
+#include "MantidSINQ/PoldiUtilities/PeakFunctionIntegrator.h"
+
+#include "MantidCurveFitting/Gaussian.h"
+#include "MantidCurveFitting/Lorentzian.h"
+
+#include "gsl/gsl_errno.h"
+
+using namespace Mantid::Poldi;
+using namespace Mantid::API;
+using namespace Mantid::CurveFitting;
+
+class PeakFunctionIntegratorTest;
+
+class TestablePeakFunctionIntegrator : public PeakFunctionIntegrator
+{
+public:
+ TestablePeakFunctionIntegrator(double requiredRelativePrecision = 1e-8) :
+ PeakFunctionIntegrator(requiredRelativePrecision)
+ {
+
+ }
+
+ friend class PeakFunctionIntegratorTest;
+};
+
+class PeakFunctionIntegratorTest : public CxxTest::TestSuite
+{
+private:
+ IPeakFunction_sptr getGaussian(double center, double fwhm, double height)
+ {
+ IPeakFunction_sptr gaussian(new Gaussian);
+ gaussian->initialize();
+ gaussian->setCentre(center);
+ gaussian->setFwhm(fwhm);
+ gaussian->setHeight(height);
+
+ return gaussian;
+ }
+
+ double getGaussianAnalyticalInfiniteIntegral(IPeakFunction_sptr gaussian)
+ {
+ return gaussian->height() * gaussian->fwhm() / (2.0 * sqrt(2.0 * log(2.0))) * sqrt(2.0 * M_PI);
+ }
+
+ IPeakFunction_sptr getLorentzian(double center, double fwhm, double height)
+ {
+ IPeakFunction_sptr lorentzian(new Lorentzian);
+ lorentzian->initialize();
+ lorentzian->setCentre(center);
+ lorentzian->setFwhm(fwhm);
+ lorentzian->setHeight(height);
+
+ return lorentzian;
+ }
+
+ double getLorentzianAnalyticalInfiniteIntegral(IPeakFunction_sptr lorentzian)
+ {
+ return lorentzian->getParameter("Amplitude");
+ }
+
+public:
+ // This pair of boilerplate methods prevent the suite being created statically
+ // This means the constructor isn't called when running other tests
+ static PeakFunctionIntegratorTest *createSuite() { return new PeakFunctionIntegratorTest(); }
+ static void destroySuite( PeakFunctionIntegratorTest *suite ) { delete suite; }
+
+ void testDefaultConstruction()
+ {
+ TestablePeakFunctionIntegrator integrator;
+
+ TS_ASSERT(integrator.m_integrationWorkspace != 0);
+ TS_ASSERT_EQUALS(integrator.m_relativePrecision, 1e-8);
+ }
+
+ void testConstruction()
+ {
+ TestablePeakFunctionIntegrator integrator(1e-10);
+
+ TS_ASSERT(integrator.m_integrationWorkspace != 0);
+ TS_ASSERT_EQUALS(integrator.m_relativePrecision, 1e-10);
+ }
+
+ void testSetRequiredRelativePrecision()
+ {
+ PeakFunctionIntegrator integrator;
+ integrator.setRequiredRelativePrecision(1e-2);
+
+ TS_ASSERT_EQUALS(integrator.requiredRelativePrecision(), 1e-2);
+ }
+
+ void testgsl_peak_wrapper()
+ {
+ IPeakFunction_sptr gaussian = getGaussian(0.0, 1.0, 2.0);
+
+ TS_ASSERT_EQUALS(gsl_peak_wrapper(0.0, &gaussian), 2.0);
+ }
+
+ void testNullPointer()
+ {
+ PeakFunctionIntegrator integrator;
+
+ TS_ASSERT_THROWS(integrator.integrateInfinity(IPeakFunction_sptr()), std::invalid_argument);
+ TS_ASSERT_THROWS(integrator.integratePositiveInfinity(IPeakFunction_sptr(), 0.0), std::invalid_argument);
+ TS_ASSERT_THROWS(integrator.integrateNegativeInfinity(IPeakFunction_sptr(), 0.0), std::invalid_argument);
+ TS_ASSERT_THROWS(integrator.integrate(IPeakFunction_sptr(), 0.0, 0.0), std::invalid_argument);
+
+ }
+
+ void testIntegrateInfinityGaussian()
+ {
+ IPeakFunction_sptr gaussian = getGaussian(0.0, 1.0, 1.0);
+
+ PeakFunctionIntegrator integrator;
+ IntegrationResult result = integrator.integrateInfinity(gaussian);
+ TS_ASSERT_EQUALS(result.errorCode, GSL_SUCCESS);
+ TS_ASSERT_DELTA(result.result, getGaussianAnalyticalInfiniteIntegral(gaussian), integrator.requiredRelativePrecision());
+ TS_ASSERT_DELTA(result.error, 0.0, integrator.requiredRelativePrecision());
+
+ integrator.setRequiredRelativePrecision(1e-14);
+ IntegrationResult otherResult = integrator.integrateInfinity(gaussian);
+ TS_ASSERT_EQUALS(otherResult.errorCode, GSL_EBADTOL);
+ TS_ASSERT_EQUALS(otherResult.result, 0.0);
+ TS_ASSERT_EQUALS(otherResult.error, 0.0);
+ }
+
+ void testIntegratePositiveInfinityGaussian()
+ {
+ IPeakFunction_sptr gaussian = getGaussian(0.0, 1.0, 1.0);
+ PeakFunctionIntegrator integrator;
+ IntegrationResult result = integrator.integratePositiveInfinity(gaussian, 0.0);
+
+ TS_ASSERT_EQUALS(result.errorCode, GSL_SUCCESS);
+ TS_ASSERT_DELTA(result.result, getGaussianAnalyticalInfiniteIntegral(gaussian) / 2.0, integrator.requiredRelativePrecision());
+ }
+
+ void testIntegrateNegativeInfinityGaussian()
+ {
+ IPeakFunction_sptr gaussian = getGaussian(0.0, 1.0, 1.0);
+ PeakFunctionIntegrator integrator;
+ IntegrationResult result = integrator.integrateNegativeInfinity(gaussian, 0.0);
+
+ TS_ASSERT_EQUALS(result.errorCode, GSL_SUCCESS);
+ TS_ASSERT_DELTA(result.result, getGaussianAnalyticalInfiniteIntegral(gaussian) / 2.0, integrator.requiredRelativePrecision());
+ }
+
+ void testIntegrateGaussian()
+ {
+ /*
+ * Normal distribution with mu = 0, sigma = 1, height = 1/sqrt(2 * pi)
+ * -integral from -1 to 1 should give approx. 0.682
+ * -integral from -2 to 2 should give approx. 0.954
+ * -integral from -3 to 3 should give approx. 0.997
+ */
+ IPeakFunction_sptr gaussian = getGaussian(0.0, 2.0 * sqrt(2.0 * log(2.0)), 1.0 / sqrt(2.0 * M_PI));
+ PeakFunctionIntegrator integrator(1e-10);
+
+ IntegrationResult rOneSigma = integrator.integrate(gaussian, -1.0, 1.0);
+ TS_ASSERT_EQUALS(rOneSigma.errorCode, GSL_SUCCESS);
+ TS_ASSERT_DELTA(rOneSigma.result, 0.682689492137086, integrator.requiredRelativePrecision());
+
+ IntegrationResult rTwoSigma = integrator.integrate(gaussian, -2.0, 2.0);
+ TS_ASSERT_EQUALS(rTwoSigma.errorCode, GSL_SUCCESS);
+ TS_ASSERT_DELTA(rTwoSigma.result, 0.954499736103642, integrator.requiredRelativePrecision());
+
+ IntegrationResult rThreeSigma = integrator.integrate(gaussian, -3.0, 3.0);
+ TS_ASSERT_EQUALS(rThreeSigma.errorCode, GSL_SUCCESS);
+ TS_ASSERT_DELTA(rThreeSigma.result, 0.997300203936740, integrator.requiredRelativePrecision());
+ }
+
+ void testIntegrateInfinityLorentzian()
+ {
+ IPeakFunction_sptr lorentzian = getLorentzian(0.0, 3.0, 8.0);
+ PeakFunctionIntegrator integrator(1e-8);
+
+ IntegrationResult result = integrator.integrateInfinity(lorentzian);
+ TS_ASSERT_EQUALS(result.errorCode, GSL_SUCCESS);
+ TS_ASSERT_DELTA(result.result, getLorentzianAnalyticalInfiniteIntegral(lorentzian), integrator.requiredRelativePrecision());
+ TS_ASSERT_LESS_THAN(result.intervals, 1000);
+ }
+
+};
+
+#endif // PEAKFUNCTIONINTEGRATORTEST_H