diff --git a/Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiFitPeaks2D.h b/Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiFitPeaks2D.h
index dc2d381822e324b24fffd5f45d4f864a9865e578..165df2ca7c0bb74e5922db6f39b58832342b1ed5 100644
--- a/Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiFitPeaks2D.h
+++ b/Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiFitPeaks2D.h
@@ -5,6 +5,10 @@
#include "MantidSINQ/DllConfig.h"
#include "MantidAPI/Algorithm.h"
#include "MantidAPI/IFunction.h"
+#include "MantidAPI/IPeakFunction.h"
+
+#include "MantidKernel/Matrix.h"
+
#include "MantidDataObjects/TableWorkspace.h"
#include "MantidSINQ/PoldiUtilities/PoldiPeakCollection.h"
#include "MantidSINQ/PoldiUtilities/PoldiTimeTransformer.h"
@@ -58,6 +62,10 @@ public:
virtual const std::string summary() const;
+ boost::shared_ptr<Kernel::DblMatrix> getLocalCovarianceMatrix(
+ const boost::shared_ptr<const Kernel::DblMatrix> &covarianceMatrix,
+ size_t parameterOffset, size_t nParams) const;
+
protected:
PoldiPeakCollection_sptr
getPeakCollection(const DataObjects::TableWorkspace_sptr &peakTable) const;
@@ -72,7 +80,7 @@ protected:
PoldiPeakCollection_sptr &to) const;
PoldiPeakCollection_sptr
- getPeakCollectionFromFunction(const API::IFunction_sptr &fitFunction) const;
+ getPeakCollectionFromFunction(const API::IFunction_sptr &fitFunction);
Poldi2DFunction_sptr getFunctionFromPeakCollection(
const PoldiPeakCollection_sptr &peakCollection) const;
void addBackgroundTerms(Poldi2DFunction_sptr poldi2DFunction) const;
diff --git a/Code/Mantid/Framework/SINQ/src/PoldiFitPeaks2D.cpp b/Code/Mantid/Framework/SINQ/src/PoldiFitPeaks2D.cpp
index 735094834f0288d4e726c83a446d37d785719773..ce552fe18c6fa3c9f00dc5b70fd85ea99ce4f156 100644
--- a/Code/Mantid/Framework/SINQ/src/PoldiFitPeaks2D.cpp
+++ b/Code/Mantid/Framework/SINQ/src/PoldiFitPeaks2D.cpp
@@ -50,9 +50,7 @@ const std::string PoldiFitPeaks2D::name() const { return "PoldiFitPeaks2D"; }
int PoldiFitPeaks2D::version() const { return 1; }
/// Algorithm's category for identification. @see Algorithm::category
-const std::string PoldiFitPeaks2D::category() const {
- return "SINQ\\Poldi";
-}
+const std::string PoldiFitPeaks2D::category() const { return "SINQ\\Poldi"; }
/// Very short algorithm summary. @see Algorith::summary
const std::string PoldiFitPeaks2D::summary() const {
@@ -102,6 +100,37 @@ void PoldiFitPeaks2D::init() {
"Table workspace with fitted peaks.");
}
+/**
+ * Extracts the covariance matrix for the supplied function
+ *
+ * This method extracts the covariance matrix for a sub-function from
+ * the global covariance matrix. If no matrix is given, a zero-matrix is
+ * returned.
+ *
+ * @param covarianceMatrix :: Global covariance matrix.
+ * @param parameterOffset :: Offset for the parameters of profileFunction.
+ * @param nParams :: Number of parameters of the local function.
+ * @return
+ */
+boost::shared_ptr<DblMatrix> PoldiFitPeaks2D::getLocalCovarianceMatrix(
+ const boost::shared_ptr<const Kernel::DblMatrix> &covarianceMatrix,
+ size_t parameterOffset, size_t nParams) const {
+ if (covarianceMatrix) {
+ boost::shared_ptr<Kernel::DblMatrix> localCov =
+ boost::make_shared<Kernel::DblMatrix>(nParams, nParams, false);
+ for (size_t j = 0; j < nParams; ++j) {
+ for (size_t k = 0; k < nParams; ++k) {
+ (*localCov)[j][k] =
+ (*covarianceMatrix)[parameterOffset + j][parameterOffset + k];
+ }
+ }
+
+ return localCov;
+ }
+
+ return boost::make_shared<Kernel::DblMatrix>(nParams, nParams, false);
+}
+
/**
* Construct a PoldiPeakCollection from a Poldi2DFunction
*
@@ -114,7 +143,7 @@ void PoldiFitPeaks2D::init() {
* @return PoldiPeakCollection containing peaks with normalized intensities
*/
PoldiPeakCollection_sptr PoldiFitPeaks2D::getPeakCollectionFromFunction(
- const IFunction_sptr &fitFunction) const {
+ const IFunction_sptr &fitFunction) {
Poldi2DFunction_sptr poldi2DFunction =
boost::dynamic_pointer_cast<Poldi2DFunction>(fitFunction);
@@ -126,6 +155,11 @@ PoldiPeakCollection_sptr PoldiFitPeaks2D::getPeakCollectionFromFunction(
PoldiPeakCollection_sptr normalizedPeaks =
boost::make_shared<PoldiPeakCollection>(PoldiPeakCollection::Integral);
+ boost::shared_ptr<const Kernel::DblMatrix> covarianceMatrix =
+ poldi2DFunction->getCovarianceMatrix();
+
+ size_t offset = 0;
+
for (size_t i = 0; i < poldi2DFunction->nFunctions(); ++i) {
boost::shared_ptr<PoldiSpectrumDomainFunction> peakFunction =
boost::dynamic_pointer_cast<PoldiSpectrumDomainFunction>(
@@ -136,36 +170,35 @@ PoldiPeakCollection_sptr PoldiFitPeaks2D::getPeakCollectionFromFunction(
boost::dynamic_pointer_cast<IPeakFunction>(
peakFunction->getProfileFunction());
+ // Get local covariance matrix
+ size_t nLocalParams = profileFunction->nParams();
+ boost::shared_ptr<Kernel::DblMatrix> localCov =
+ getLocalCovarianceMatrix(covarianceMatrix, offset, nLocalParams);
+ profileFunction->setCovarianceMatrix(localCov);
+
+ // Increment offset for next function
+ offset += nLocalParams;
+
+ IAlgorithm_sptr errorAlg = createChildAlgorithm("EstimatePeakErrors");
+ errorAlg->setProperty(
+ "Function", boost::dynamic_pointer_cast<IFunction>(profileFunction));
+ errorAlg->setPropertyValue("OutputWorkspace", "Errors");
+ errorAlg->execute();
+
double centre = profileFunction->centre();
double height = profileFunction->height();
+ double fwhmValue = profileFunction->fwhm();
- size_t dIndex = 0;
- size_t iIndex = 0;
- size_t fIndex = 0;
-
- for (size_t j = 0; j < profileFunction->nParams(); ++j) {
- if (profileFunction->getParameter(j) == centre) {
- dIndex = j;
- } else if (profileFunction->getParameter(j) == height) {
- iIndex = j;
- } else {
- fIndex = j;
- }
- }
-
- // size_t dIndex = peakFunction->parameterIndex("Centre");
- UncertainValue d(peakFunction->getParameter(dIndex),
- peakFunction->getError(dIndex));
+ ITableWorkspace_sptr errorTable =
+ errorAlg->getProperty("OutputWorkspace");
- // size_t iIndex = peakFunction->parameterIndex("Area");
- UncertainValue intensity(peakFunction->getParameter(iIndex),
- peakFunction->getError(iIndex));
+ double centreError = errorTable->cell<double>(0, 2);
+ double heightError = errorTable->cell<double>(1, 2);
+ double fwhmError = errorTable->cell<double>(2, 2);
- // size_t fIndex = peakFunction->parameterIndex("Sigma");
- double fwhmValue = profileFunction->fwhm();
- UncertainValue fwhm(fwhmValue, fwhmValue /
- peakFunction->getParameter(fIndex) *
- peakFunction->getError(fIndex));
+ UncertainValue d(centre, centreError);
+ UncertainValue intensity(height, heightError);
+ UncertainValue fwhm(fwhmValue, fwhmError);
PoldiPeak_sptr peak =
PoldiPeak::create(MillerIndices(), d, intensity, UncertainValue(1.0));
diff --git a/Code/Mantid/Framework/SINQ/test/PoldiFitPeaks2DTest.h b/Code/Mantid/Framework/SINQ/test/PoldiFitPeaks2DTest.h
index 9c18d599c85b0f0833c8ea46ca5a989824390667..1ac7a69edaaf3814ff1e7001da18d7d75adc1249 100644
--- a/Code/Mantid/Framework/SINQ/test/PoldiFitPeaks2DTest.h
+++ b/Code/Mantid/Framework/SINQ/test/PoldiFitPeaks2DTest.h
@@ -6,6 +6,8 @@
#include "MantidAPI/FrameworkManager.h"
#include "MantidAPI/AlgorithmManager.h"
+#include "MantidKernel/Matrix.h"
+
#include "MantidSINQ/PoldiFitPeaks2D.h"
#include "MantidSINQ/PoldiUtilities/PoldiSpectrumDomainFunction.h"
#include "MantidSINQ/PoldiUtilities/PoldiMockInstrumentHelpers.h"
@@ -14,6 +16,7 @@
using namespace Mantid::Poldi;
using namespace Mantid::API;
using namespace Mantid::DataObjects;
+using namespace Mantid::Kernel;
class PoldiFitPeaks2DTest : public CxxTest::TestSuite
{
@@ -234,6 +237,18 @@ public:
PoldiPeakCollection_sptr peaks = PoldiPeakCollectionHelpers::createPoldiPeakCollectionNormalized();
IFunction_sptr poldi2DFunction = spectrumCalculator.getFunctionFromPeakCollection(peaks);
+ size_t nParams = poldi2DFunction->nParams();
+
+ // Make a matrix with diagonal elements = 0.05 and set as covariance matrix
+ boost::shared_ptr<DblMatrix> matrix = boost::make_shared<DblMatrix>(nParams, nParams, true);
+ matrix->operator *=(0.05);
+ poldi2DFunction->setCovarianceMatrix(matrix);
+
+ // Also set errors for old behavior
+ for(size_t i = 0; i < nParams; ++i) {
+ poldi2DFunction->setError(i, sqrt(0.05));
+ }
+
PoldiPeakCollection_sptr peaksFromFunction = spectrumCalculator.getPeakCollectionFromFunction(poldi2DFunction);
TS_ASSERT_EQUALS(peaksFromFunction->peakCount(), peaks->peakCount());
@@ -241,8 +256,11 @@ public:
PoldiPeak_sptr functionPeak = peaksFromFunction->peak(i);
PoldiPeak_sptr referencePeak = peaks->peak(i);
- TS_ASSERT_EQUALS(functionPeak->d(), referencePeak->d());
- TS_ASSERT_EQUALS(functionPeak->fwhm(), referencePeak->fwhm());
+ TS_ASSERT_EQUALS(functionPeak->d().value(), referencePeak->d().value());
+ TS_ASSERT_EQUALS(functionPeak->fwhm().value(), referencePeak->fwhm().value());
+ TS_ASSERT_DELTA(functionPeak->d().error(), sqrt(0.05), 1e-6);
+ TS_ASSERT_DELTA(functionPeak->fwhm(PoldiPeak::AbsoluteD).error(), sqrt(0.05) * (2.0 * sqrt(2.0 * log(2.0))), 1e-6);
+
}
}