VesuvioResolution.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/Functions/VesuvioResolution.h"
#include "MantidAPI/FunctionFactory.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/SpectrumInfo.h"
#include "MantidCurveFitting/Algorithms/ConvertToYSpace.h"
#include "MantidGeometry/Instrument.h"
#include "MantidKernel/PhysicalConstants.h"
 
#include <gsl/gsl_poly.h>
 
namespace Mantid::CurveFitting::Functions {
 
using namespace Mantid;
using namespace Mantid::CurveFitting;
using namespace Mantid::CurveFitting::Algorithms;
namespace {
const char *MASS_NAME = "Mass";
 
const double STDDEV_TO_HWHM = std::sqrt(std::log(4.0));
} // namespace
 
// Register into factory
DECLARE_FUNCTION(VesuvioResolution)
 
 
ResolutionParams VesuvioResolution::getResolutionParameters(const API::MatrixWorkspace_const_sptr &ws,
                                                            const size_t index) {
  const auto &spectrumInfo = ws->spectrumInfo();
  if (!spectrumInfo.hasDetectors(index))
    throw std::invalid_argument("VesuvioResolution - Workspace has no detector "
                                "attached to histogram at index " +
                                std::to_string(index));
 
  ResolutionParams respar;
  const auto &pmap = ws->constInstrumentParameters();
  const auto &det = spectrumInfo.detector(index);
  respar.dl1 = ConvertToYSpace::getComponentParameter(det, pmap, "sigma_l1");
  respar.dl2 = ConvertToYSpace::getComponentParameter(det, pmap, "sigma_l2");
  respar.dtof = ConvertToYSpace::getComponentParameter(det, pmap, "sigma_tof");
  respar.dthe = ConvertToYSpace::getComponentParameter(det, pmap, "sigma_theta"); // radians
  respar.dEnLorentz = ConvertToYSpace::getComponentParameter(det, pmap, "hwhm_lorentz");
  respar.dEnGauss = ConvertToYSpace::getComponentParameter(det, pmap, "sigma_gauss");
  return respar;
}
 
VesuvioResolution::VesuvioResolution()
    : API::ParamFunction(), API::IFunction1D(), m_log("VesuvioResolution"), m_wsIndex(0), m_mass(0.0), m_voigt(),
      m_resolutionSigma(0.0), m_lorentzFWHM(0.0) {}
 
std::string VesuvioResolution::name() const { return "VesuvioResolution"; }
 
/*
 * Creates the internal caches
 */
void VesuvioResolution::setUpForFit() {
  // Voigt
 
  using namespace Mantid::API;
  m_voigt = std::dynamic_pointer_cast<IPeakFunction>(FunctionFactory::Instance().createFunction("Voigt"));
}
 
void VesuvioResolution::setMatrixWorkspace(std::shared_ptr<const API::MatrixWorkspace> workspace, size_t wsIndex,
                                           double startX, double endX) {
  UNUSED_ARG(startX);
  UNUSED_ARG(endX);
 
  m_wsIndex = wsIndex;
  DetectorParams detpar = ConvertToYSpace::getDetectorParameters(workspace, m_wsIndex);
  ResolutionParams respar = getResolutionParameters(workspace, m_wsIndex);
  this->cacheResolutionComponents(detpar, respar);
}
 
void VesuvioResolution::cacheResolutionComponents(const Algorithms::DetectorParams &detpar,
                                                  const ResolutionParams &respar) {
  // geometry
  double theta = detpar.theta; // cache for frequent access
  double hwhmLorentzE = respar.dEnLorentz;
  double hwhmGaussE = STDDEV_TO_HWHM * respar.dEnGauss;
 
  // ------ Fixed coefficients related to resolution & Y-space transforms
  // ------------------
  const double mn = PhysicalConstants::NeutronMassAMU;
  const double mevToK = PhysicalConstants::E_mev_toNeutronWavenumberSq;
 
  const double k1 = std::sqrt(detpar.efixed / mevToK);
  const double l2l1 = detpar.l2 / detpar.l1;
 
  // Resolution dependence
 
  // Find K0/K1 at y=0 by taking the largest root of (M-1)s^2 + 2cos(theta)s -
  // (M+1) = 0
  // Quadratic if M != 1 but simple linear if it does
  double k0k1(0.0);
  if ((m_mass - 1.0) > DBL_EPSILON) {
    double x0(0.0), x1(0.0);
    gsl_poly_solve_quadratic(m_mass - 1.0, 2.0 * std::cos(theta), -(m_mass + 1.0), &x0, &x1);
    k0k1 = std::max(x0, x1); // K0/K1 at y=0
  } else {
    // solution is simply s = 1/cos(theta)
    k0k1 = 1.0 / std::cos(theta);
  }
  double qy0(0.0), wgauss(0.0);
 
  if (m_mass > 1.0) {
    qy0 = std::sqrt(k1 * k1 * m_mass * (k0k1 * k0k1 - 1));
    double k0k1p3 = std::pow(k0k1, 3);
    double r1 = -(1.0 + l2l1 * k0k1p3);
    double r2 = 1.0 - l2l1 * k0k1p3 + l2l1 * std::pow(k0k1, 2) * std::cos(theta) - k0k1 * std::cos(theta);
 
    double factor = (0.2413 / qy0) * ((m_mass / mn) * r1 - r2);
    m_lorentzFWHM = std::abs(factor * hwhmLorentzE * 2);
    wgauss = std::abs(factor * hwhmGaussE * 2);
  } else {
    qy0 = k1 * std::tan(theta);
    double factor = (0.2413 * 2.0 / k1) * std::abs((std::cos(theta) + l2l1) / std::sin(theta));
    m_lorentzFWHM = hwhmLorentzE * factor;
    wgauss = hwhmGaussE * factor;
  }
 
  double k0y0 = k1 * k0k1; // k0_y0 =  k0 value at y=0
  double wtheta = 2.0 * STDDEV_TO_HWHM * std::abs(k0y0 * k1 * std::sin(theta) / qy0) * respar.dthe;
  double common = (m_mass / mn) - 1 + k1 * std::cos(theta) / k0y0;
  double wl1 = 2.0 * STDDEV_TO_HWHM * std::abs((std::pow(k0y0, 2) / (qy0 * detpar.l1)) * common) * respar.dl1;
  double wl2 = 2.0 * STDDEV_TO_HWHM * std::abs((std::pow(k0y0, 3) / (k1 * qy0 * detpar.l1)) * common) * respar.dl2;
 
  m_resolutionSigma = std::sqrt(std::pow(wgauss, 2) + std::pow(wtheta, 2) + std::pow(wl1, 2) + std::pow(wl2, 2));
 
  m_log.information() << "--------------------- Mass=" << m_mass << " -----------------------\n";
  m_log.information() << "w_l1 (FWHM)=" << wl2 << '\n';
  m_log.information() << "w_l0 (FWHM)=" << wl1 << '\n';
  m_log.information() << "w_theta (FWHM)=" << wtheta << '\n';
  m_log.information() << "w_foil_lorentz (FWHM)=" << m_lorentzFWHM << '\n';
  m_log.information() << "w_foil_gauss (FWHM)=" << wgauss << '\n';
}
 
void VesuvioResolution::function1D(double *out, const double *xValues, const size_t nData) const {
  std::vector<double> outVec(static_cast<int>(nData), 0);
  const std::vector<double> xValuesVec(xValues, xValues + nData);
  voigtApprox(outVec, xValuesVec, 0, 1, m_lorentzFWHM, m_resolutionSigma);
  std::copy(outVec.begin(), outVec.end(), out);
}
 
void VesuvioResolution::declareAttributes() { declareAttribute(MASS_NAME, IFunction::Attribute(m_mass)); }
 
void VesuvioResolution::setAttribute(const std::string &name, const Attribute &value) {
  IFunction::setAttribute(name, value); // Make sure the base-class stores it
  if (name == MASS_NAME)
    m_mass = value.asDouble();
}
 
void VesuvioResolution::voigtApprox(std::vector<double> &voigt, const std::vector<double> &xValues,
                                    const double lorentzPos, const double lorentzAmp) const {
  voigtApprox(voigt, xValues, lorentzPos, lorentzAmp, m_lorentzFWHM, m_resolutionSigma);
}
 
void VesuvioResolution::voigtApprox(std::vector<double> &voigt, const std::vector<double> &xValues,
                                    const double lorentzPos, const double lorentzAmp, const double lorentzWidth,
                                    const double gaussWidth) const {
  m_voigt->setParameter(0, lorentzAmp);
  m_voigt->setParameter(1, lorentzPos);
  m_voigt->setParameter(2, lorentzWidth);
  m_voigt->setParameter(3, gaussWidth);
  assert(voigt.size() == xValues.size());
  m_voigt->functionLocal(voigt.data(), xValues.data(), xValues.size());
 
  // Normalize so that integral of V=lorentzAmp
  const double norm = 1.0 / (0.5 * M_PI * lorentzWidth);
  using std::placeholders::_1;
  std::transform(voigt.begin(), voigt.end(), voigt.begin(), std::bind(std::multiplies<double>(), _1, norm));
}
 
} // namespace Mantid::CurveFitting::Functions