EstimatePeakErrors.cpp

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidCurveFitting/Algorithms/EstimatePeakErrors.h"
#include "MantidCurveFitting/EigenMatrix.h"
#include "MantidCurveFitting/Functions/PeakParameterFunction.h"
 
#include "MantidAPI/CompositeFunction.h"
#include "MantidAPI/FunctionProperty.h"
#include "MantidAPI/IPeakFunction.h"
#include "MantidAPI/ITableWorkspace.h"
#include "MantidAPI/TableRow.h"
#include "MantidAPI/WorkspaceFactory.h"
 
namespace Mantid::CurveFitting::Algorithms {
using namespace API;
using namespace Kernel;
using namespace std;
 
// Subscription
DECLARE_ALGORITHM(EstimatePeakErrors)
 
//--------------------------------------------------------------------------------------------------------
// Public members
//--------------------------------------------------------------------------------------------------------
 
 
EstimatePeakErrors::EstimatePeakErrors() : Algorithm() {}
 
const std::string EstimatePeakErrors::summary() const {
  return "Calculates error estimates for peak parameters: "
         "centre, height, FWHM and intensity.";
}
 
const std::string EstimatePeakErrors::name() const { return "EstimatePeakErrors"; }
 
int EstimatePeakErrors::version() const { return 1; }
 
const std::string EstimatePeakErrors::category() const { return "Optimization"; }
 
//--------------------------------------------------------------------------------------------------------
// Private members
//--------------------------------------------------------------------------------------------------------
namespace {
 
EigenMatrix makeJacobian(IPeakFunction &peak, double &centre, double &height, double &fwhm, double &intensity) {
  EigenMatrix jacobian(4, peak.nParams());
  centre = peak.centre();
  height = peak.height();
  fwhm = peak.fwhm();
  intensity = peak.intensity();
  for (size_t ip = 0; ip < peak.nParams(); ++ip) {
    double p = peak.getParameter(ip);
    double dp = 1e-9;
    if (p != 0.0) {
      dp *= p;
    }
    peak.setParameter(ip, p + dp);
    jacobian.set(0, ip, (peak.centre() - centre) / dp);
    jacobian.set(1, ip, (peak.height() - height) / dp);
    jacobian.set(2, ip, (peak.fwhm() - fwhm) / dp);
    jacobian.set(3, ip, (peak.intensity() - intensity) / dp);
    peak.setParameter(ip, p);
  }
  return jacobian;
}
 
void calculatePeakValues(IPeakFunction &peak, ITableWorkspace &results, const EigenMatrix &covariance,
                         const std::string &prefix) {
  double centre, height, fwhm, intensity;
  EigenMatrix J = makeJacobian(peak, centre, height, fwhm, intensity);
  EigenMatrix JCJ = J * covariance * J.tr();
 
  TableRow row = results.appendRow();
  row << prefix + "Centre" << centre << sqrt(JCJ.get(0, 0));
  row = results.appendRow();
  row << prefix + "Height" << height << sqrt(JCJ.get(1, 1));
  row = results.appendRow();
  row << prefix + "FWHM" << fwhm << sqrt(JCJ.get(2, 2));
  row = results.appendRow();
  row << prefix + "Intensity" << intensity << sqrt(JCJ.get(3, 3));
}
} // namespace
 
void EstimatePeakErrors::init() {
 
  declareProperty(std::make_unique<FunctionProperty>("Function", Direction::InOut),
                  "Fitting function containing peaks. Must have a covariance "
                  "matrix attached.");
 
  declareProperty(
      std::make_unique<API::WorkspaceProperty<API::ITableWorkspace>>("OutputWorkspace", "", Kernel::Direction::Output),
      "The name of the TableWorkspace with the output values and errors.");
}
 
void EstimatePeakErrors::exec() {
 
  IFunction_sptr function = getProperty("Function");
 
  ITableWorkspace_sptr results = WorkspaceFactory::Instance().createTable("TableWorkspace");
  results->addColumn("str", "Parameter");
  results->addColumn("double", "Value");
  results->addColumn("double", "Error");
 
  auto matrix = function->getCovarianceMatrix();
  if (!matrix) {
    g_log.warning() << "Function doesn't have covariance matrix.\n";
    setProperty("OutputWorkspace", results);
    return;
  }
 
  auto *peak = dynamic_cast<IPeakFunction *>(function.get());
 
  if (peak) {
    EigenMatrix covariance(*matrix);
    calculatePeakValues(*peak, *results, covariance, "");
  } else {
    auto *cf = dynamic_cast<CompositeFunction *>(function.get());
    if (cf) {
      size_t ip = 0;
      for (size_t i = 0; i < cf->nFunctions(); ++i) {
        IFunction *fun = cf->getFunction(i).get();
        size_t np = fun->nParams();
        peak = dynamic_cast<IPeakFunction *>(fun);
        if (peak) {
          std::string prefix = "f" + std::to_string(i) + ".";
          EigenMatrix covariance(*matrix, ip, ip, np, np);
          calculatePeakValues(*peak, *results, covariance, prefix);
        }
        ip += np;
      }
    } else {
      g_log.warning() << "Function has no peaks.\n";
    }
  }
 
  setProperty("OutputWorkspace", results);
}
 
} // namespace Mantid::CurveFitting::Algorithms