GaussianComptonProfile.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 +
//------------------------------------------------------------------------------------------------
// Includes
//------------------------------------------------------------------------------------------------
#include "MantidCurveFitting/Functions/GaussianComptonProfile.h"
#include "MantidAPI/FunctionFactory.h"
 
#include <cmath>
 
namespace Mantid::CurveFitting::Functions {
 
using namespace CurveFitting;
DECLARE_FUNCTION(GaussianComptonProfile)
 
const char *WIDTH_PARAM = "Width";
const char *AMP_PARAM = "Intensity";
 
const double STDDEV_TO_HWHM = std::sqrt(std::log(4.0));
 
GaussianComptonProfile::GaussianComptonProfile() : ComptonProfile() {}
 
std::string GaussianComptonProfile::name() const { return "GaussianComptonProfile"; }
 
void GaussianComptonProfile::declareParameters() {
  ComptonProfile::declareParameters();
  declareParameter(WIDTH_PARAM, 1.0, "Gaussian width parameter");
  declareParameter(AMP_PARAM, 1.0, "Gaussian intensity parameter");
}
 
std::vector<size_t> GaussianComptonProfile::intensityParameterIndices() const {
  return std::vector<size_t>(1, this->parameterIndex(AMP_PARAM));
}
 
size_t GaussianComptonProfile::fillConstraintMatrix(Kernel::DblMatrix &cmatrix, const size_t start,
                                                    const HistogramData::HistogramE &errors) const {
  std::vector<double> result(ySpace().size());
  const double amplitude = 1.0;
  this->massProfile(result.data(), ySpace().size(), amplitude);
  std::transform(result.begin(), result.end(), errors.begin(), result.begin(), std::divides<double>());
  cmatrix.setColumn(start, result);
  return 1;
}
 
void GaussianComptonProfile::massProfile(double *result, const size_t nData) const {
  const double amplitude(getParameter(AMP_PARAM));
  this->massProfile(result, nData, amplitude);
}
 
void GaussianComptonProfile::massProfile(double *result, const size_t nData, const double amplitude) const {
  double lorentzPos(0.0), gaussWidth(getParameter(WIDTH_PARAM));
  double gaussFWHM =
      std::sqrt(std::pow(m_resolutionFunction->resolutionFWHM(), 2) + std::pow(2.0 * STDDEV_TO_HWHM * gaussWidth, 2));
 
  const auto &yspace = ySpace();
  // Gaussian already folded into Voigt
  std::vector<double> voigt(yspace.size()), voigtDiffResult(yspace.size());
  m_resolutionFunction->voigtApprox(voigt, yspace, lorentzPos, amplitude, m_resolutionFunction->lorentzFWHM(),
                                    gaussFWHM);
  voigtApproxDiff(voigtDiffResult, yspace, lorentzPos, amplitude, m_resolutionFunction->lorentzFWHM(), gaussFWHM);
 
  const auto &modq = modQ();
  const auto &ei = e0();
  if (modq.empty() || ei.empty()) {
    throw std::runtime_error("The Q values or e0 values have not been set");
  }
  // Include e_i^0.1*mass/q pre-factor
  for (size_t j = 0; j < nData; ++j) {
    const double q = modq[j];
    const double prefactor = mass() * std::pow(ei[j], 0.1) / q;
    result[j] = prefactor * (voigt[j] - std::pow(gaussWidth, 4.0) * voigtDiffResult[j] / (3.0 * q));
  }
}
 
} // namespace Mantid::CurveFitting::Functions