Merge pull request #240 from mantidproject/11104_add_intensity_to_ipe…

…akfunction

Add intensity/setIntensity methods to IPeakFunction

Code/Mantid/Framework/API/CMakeLists.txt

@@ -93,6 +93,7 @@ set ( SRC_FILES
 	src/ParamFunction.cpp
 	src/ParameterReference.cpp
 	src/ParameterTie.cpp
+        src/PeakFunctionIntegrator.cpp
 	src/PeakTransform.cpp
 	src/PeakTransformHKL.cpp
 	src/PeakTransformQLab.cpp
@@ -252,6 +253,7 @@ set ( INC_FILES
 	inc/MantidAPI/ParamFunction.h
 	inc/MantidAPI/ParameterReference.h
 	inc/MantidAPI/ParameterTie.h
+        inc/MantidAPI/PeakFunctionIntegrator.h
 	inc/MantidAPI/PeakTransform.h
 	inc/MantidAPI/PeakTransformFactory.h
 	inc/MantidAPI/PeakTransformHKL.h
@@ -348,6 +350,7 @@ set ( TEST_FILES
 	ParamFunctionAttributeHolderTest.h
 	ParameterReferenceTest.h
 	ParameterTieTest.h
+        PeakFunctionIntegratorTest.h
 	PeakTransformHKLTest.h
 	PeakTransformQLabTest.h
 	PeakTransformQSampleTest.h

Code/Mantid/Framework/API/inc/MantidAPI/IPeakFunction.h

@@ -49,6 +49,12 @@ class MANTID_API_DLL IPeakFunction : public IFunctionWithLocation {
   /// Sets the parameters such that FWHM = w
   virtual void setFwhm(const double w) = 0;
 
+  /// Returns the integral intensity of the peak
+  virtual double intensity() const;
+
+  /// Sets the integral intensity of the peak
+  virtual void setIntensity(const double newIntensity);
+
   /// General implementation of the method for all peaks.
   void function1D(double *out, const double *xValues, const size_t nData) const;
   /// General implementation of the method for all peaks.

Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiUtilities/PeakFunctionIntegrator.h → Code/Mantid/Framework/API/inc/MantidAPI/PeakFunctionIntegrator.h

@@ -1,22 +1,25 @@
 #ifndef PEAKFUNCTIONINTEGRATOR_H
 #define PEAKFUNCTIONINTEGRATOR_H
 
-#include "MantidSINQ/DllConfig.h"
+#include "MantidAPI/DllConfig.h"
 #include "MantidAPI/IPeakFunction.h"
 #include "gsl/gsl_integration.h"
 
 namespace Mantid {
-namespace Poldi {
+namespace API {
 
 /** PeakFunctionIntegrator :
  *
   General integration of peaks (in the form of IPeakFunction) by wrapping the
   corresponding GSL-functions. Integration with infinity limits is supported.
 
+  PeakFunctionIntegrator allocates a GSL integration workspace on construction
+  and frees the memory when it's destroyed.
+
     @author Michael Wedel, Paul Scherrer Institut - SINQ
     @date 24/04/2014
 
-    Copyright © 2014 PSI-MSS
+    Copyright © 2014,2015 PSI-MSS
 
   This file is part of Mantid.
 
@@ -37,7 +40,7 @@ namespace Poldi {
   Code Documentation is available at: <http://doxygen.mantidproject.org>
 */
 
-struct MANTID_SINQ_DLL IntegrationResult {
+struct MANTID_API_DLL IntegrationResult {
   double result;
   double error;
   size_t intervals;
@@ -46,7 +49,7 @@ struct MANTID_SINQ_DLL IntegrationResult {
   bool success;
 };
 
-class MANTID_SINQ_DLL PeakFunctionIntegrator {
+class MANTID_API_DLL PeakFunctionIntegrator {
 public:
   PeakFunctionIntegrator(double requiredRelativePrecision = 1e-8);
   virtual ~PeakFunctionIntegrator();
@@ -55,27 +58,26 @@ class MANTID_SINQ_DLL PeakFunctionIntegrator {
   double requiredRelativePrecision() const;
 
   IntegrationResult
-  integrateInfinity(API::IPeakFunction_const_sptr peakFunction) const;
+  integrateInfinity(const IPeakFunction &peakFunction) const;
   IntegrationResult
-  integratePositiveInfinity(API::IPeakFunction_const_sptr peakFunction,
+  integratePositiveInfinity(const IPeakFunction &peakFunction,
                             double lowerLimit) const;
   IntegrationResult
-  integrateNegativeInfinity(API::IPeakFunction_const_sptr peakFunction,
+  integrateNegativeInfinity(const IPeakFunction &peakFunction,
                             double upperLimit) const;
 
-  IntegrationResult integrate(API::IPeakFunction_const_sptr peakFunction,
+  IntegrationResult integrate(const IPeakFunction &peakFunction,
                               double lowerLimit, double upperLimit) const;
 
 protected:
-  gsl_function getGSLFunction(API::IPeakFunction_const_sptr peakFunction) const;
-  void throwIfInvalid(API::IPeakFunction_const_sptr peakFunction) const;
+  gsl_function getGSLFunction(const IPeakFunction &peakFunction) const;
 
   gsl_integration_workspace *m_integrationWorkspace;
 
   double m_relativePrecision;
 };
 
-double MANTID_SINQ_DLL gsl_peak_wrapper(double x, void *parameters);
+double MANTID_API_DLL gsl_peak_wrapper(double x, void *parameters);
 }
 }
 

Code/Mantid/Framework/API/src/IPeakFunction.cpp

@@ -3,6 +3,7 @@
 //----------------------------------------------------------------------
 #include "MantidAPI/IPeakFunction.h"
 #include "MantidAPI/Jacobian.h"
+#include "MantidAPI/PeakFunctionIntegrator.h"
 #include "MantidKernel/Exception.h"
 #include "MantidKernel/ConfigService.h"
 
@@ -130,5 +131,43 @@ void IPeakFunction::setPeakRadius(const int &r) {
   }
 }
 
+/// Returns the integral intensity of the peak function, using the peak radius
+/// to determine integration borders.
+double IPeakFunction::intensity() const {
+  double x0 = centre();
+  double dx = fabs(s_peakRadius * fwhm());
+
+  PeakFunctionIntegrator integrator;
+  IntegrationResult result = integrator.integrate(*this, x0 - dx, x0 + dx);
+
+  if (!result.success) {
+    return 0.0;
+  }
+
+  return result.result;
+}
+
+/// Sets the integral intensity of the peak by adjusting the height.
+void IPeakFunction::setIntensity(const double newIntensity) {
+  double currentHeight = height();
+  double currentIntensity = intensity();
+
+  if (currentIntensity == 0.0) {
+    // Try to set a different height first.
+    setHeight(2.0);
+
+    currentHeight = height();
+    currentIntensity = intensity();
+
+    // If the current intensity is still 0, there's nothing left to do.
+    if (currentIntensity == 0.0) {
+      throw std::invalid_argument(
+          "Cannot set new intensity, not enough information available.");
+    }
+  }
+
+  setHeight(newIntensity / currentIntensity * currentHeight);
+}
+
 } // namespace API
 } // namespace Mantid

Code/Mantid/Framework/SINQ/src/PoldiUtilities/PeakFunctionIntegrator.cpp → Code/Mantid/Framework/API/src/PeakFunctionIntegrator.cpp

@@ -1,14 +1,12 @@
-#include "MantidSINQ/PoldiUtilities/PeakFunctionIntegrator.h"
+#include "MantidAPI/PeakFunctionIntegrator.h"
 
 #include "MantidAPI/FunctionDomain1D.h"
 #include "gsl/gsl_errno.h"
 #include <iostream>
 #include <iomanip>
 
 namespace Mantid {
-namespace Poldi {
-
-using namespace API;
+namespace API {
 
 /** Constructor with required relative precision argument. The default is 1e-8.
  *  See also PeakFunctionIntegrator::setRequiredRelativePrecision.
@@ -31,8 +29,7 @@ PeakFunctionIntegrator::~PeakFunctionIntegrator() {
 }
 
 /** This method sets the desired numerical relative precision that's passed on
- *to the
- *  GSL integration-routines.
+ *  to the GSL integration-routines.
  *
  *  @param newPrecision :: Desired relative precision for integrations.
  */
@@ -47,19 +44,14 @@ double PeakFunctionIntegrator::requiredRelativePrecision() const {
 }
 
 /** 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).
+ *  gsl_integration_qagi is used for this. 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);
-
+    const IPeakFunction &peakFunction) const {
   IntegrationResult result;
 
   gsl_function f = getGSLFunction(peakFunction);
@@ -74,16 +66,12 @@ IntegrationResult PeakFunctionIntegrator::integrateInfinity(
 }
 
 /** 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.
+ *  gsl_integration_qagiu is used for this.
  *
  *  @param peakFunction :: Peak function to integrate.
  */
 IntegrationResult PeakFunctionIntegrator::integratePositiveInfinity(
-    IPeakFunction_const_sptr peakFunction, double lowerLimit) const {
-  throwIfInvalid(peakFunction);
-
+    const IPeakFunction &peakFunction, double lowerLimit) const {
   IntegrationResult result;
 
   gsl_function f = getGSLFunction(peakFunction);
@@ -98,16 +86,12 @@ IntegrationResult PeakFunctionIntegrator::integratePositiveInfinity(
 }
 
 /** 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.
+ *  gsl_integration_qagil is used for this.
  *
  *  @param peakFunction :: Peak function to integrate.
  */
 IntegrationResult PeakFunctionIntegrator::integrateNegativeInfinity(
-    IPeakFunction_const_sptr peakFunction, double upperLimit) const {
-  throwIfInvalid(peakFunction);
-
+    const IPeakFunction &peakFunction, double upperLimit) const {
   IntegrationResult result;
 
   gsl_function f = getGSLFunction(peakFunction);
@@ -122,17 +106,13 @@ IntegrationResult PeakFunctionIntegrator::integrateNegativeInfinity(
 }
 
 /** 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.
+ *  gsl_integration_qags is used for this.
  *
  *  @param peakFunction :: Peak function to integrate.
  */
 IntegrationResult
-PeakFunctionIntegrator::integrate(IPeakFunction_const_sptr peakFunction,
+PeakFunctionIntegrator::integrate(const IPeakFunction &peakFunction,
                                   double lowerLimit, double upperLimit) const {
-  throwIfInvalid(peakFunction);
-
   IntegrationResult result;
 
   gsl_function f = getGSLFunction(peakFunction);
@@ -151,24 +131,22 @@ PeakFunctionIntegrator::integrate(IPeakFunction_const_sptr peakFunction,
  *  @param peakFunction :: Peak function to wrap.
  */
 gsl_function PeakFunctionIntegrator::getGSLFunction(
-    IPeakFunction_const_sptr peakFunction) const {
+    const IPeakFunction &peakFunction) const {
   gsl_function f;
-  f.function = &Mantid::Poldi::gsl_peak_wrapper;
-  f.params = &peakFunction;
+  f.function = &Mantid::API::gsl_peak_wrapper;
+  f.params =
+      reinterpret_cast<void *>(&const_cast<IPeakFunction &>(peakFunction));
 
   return f;
 }
 
-void PeakFunctionIntegrator::throwIfInvalid(
-    IPeakFunction_const_sptr peakFunction) const {
+double gsl_peak_wrapper(double x, void *parameters) {
+  IPeakFunction *peakFunction = reinterpret_cast<IPeakFunction *>(parameters);
+
   if (!peakFunction) {
-    throw std::invalid_argument("Can not integrate NULL-function.");
+    throw std::runtime_error(
+        "Cannot process NULL-pointer in gsl_peak_wrapper.");
   }
-}
-
-double gsl_peak_wrapper(double x, void *parameters) {
-  IPeakFunction_const_sptr peakFunction =
-      *(IPeakFunction_const_sptr *)parameters;
 
   double y;