InelasticDiffRotDiscreteCircle.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/InelasticDiffRotDiscreteCircle.h"
#include "MantidCurveFitting/Constraints/BoundaryConstraint.h"
 
#include "MantidAPI/FunctionFactory.h"
#include "MantidAPI/IFunction.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/SpectrumInfo.h"
#include "MantidGeometry/IDetector.h"
#include "MantidGeometry/Instrument/DetectorInfo.h"
#include "MantidKernel/Exception.h"
#include "MantidKernel/UnitConversion.h"
 
#include "MantidTypes/SpectrumDefinition.h"
 
#include <cmath>
#include <limits>
#include <sstream>
 
using BConstraint = Mantid::CurveFitting::Constraints::BoundaryConstraint;
 
namespace {
Mantid::Kernel::Logger g_log("InelasticDiffDiffRotDiscreteCircle");
}
 
namespace Mantid::CurveFitting::Functions {
 
DECLARE_FUNCTION(InelasticDiffRotDiscreteCircle)
 
 
InelasticDiffRotDiscreteCircle::InelasticDiffRotDiscreteCircle() : m_hbar(0.658211626) {
  this->declareParameter("Intensity", 1.0, "scaling factor [no units]");
  this->declareParameter("Radius", 1.0, "Circle radius [Angstroms]");
  this->declareParameter("Decay", 1.0,
                         "Inverse of transition rate, in nanoseconds "
                         "if energy in micro-ev, or picoseconds if "
                         "energy in mili-eV");
  this->declareParameter("Shift", 0.0, "Shift in the centre of the peak");
 
  this->declareAttribute("N", API::IFunction::Attribute(3));
  declareAttributes();
}
 
void InelasticDiffRotDiscreteCircle::init() {
  // Ensure positive values for Intensity, Radius, and decay
  auto IntensityConstraint =
      std::make_unique<BConstraint>(this, "Intensity", std::numeric_limits<double>::epsilon(), true);
  this->addConstraint(std::move(IntensityConstraint));
 
  auto RadiusConstraint = std::make_unique<BConstraint>(this, "Radius", std::numeric_limits<double>::epsilon(), true);
  this->addConstraint(std::move(RadiusConstraint));
 
  auto DecayConstraint = std::make_unique<BConstraint>(this, "Decay", std::numeric_limits<double>::epsilon(), true);
  this->addConstraint(std::move(DecayConstraint));
}
 
void InelasticDiffRotDiscreteCircle::function1D(double *out, const double *xValues, const size_t nData) const {
 
  auto I = this->getParameter("Intensity");
  auto R = this->getParameter("Radius");
  auto rate = m_hbar / this->getParameter("Decay"); // micro-eV or mili-eV
  auto N = this->getAttribute("N").asInt();
 
  // Retrieve Q-value from the appropriate attribute
  double Q;
  if (this->getAttribute("Q").asDouble() == EMPTY_DBL()) {
    if (m_qValueCache.empty()) {
      throw std::runtime_error("No Q attribute provided and cannot retrieve from workspace.");
    }
    const int specIdx = this->getAttribute("WorkspaceIndex").asInt();
    Q = m_qValueCache[specIdx];
 
    g_log.debug() << "Get Q value for workspace index " << specIdx << ": " << Q << '\n';
  } else {
    Q = getAttribute("Q").asDouble();
 
    g_log.debug() << "Using Q attribute: " << Q << '\n';
  }
 
  std::vector<double> sph(N);
  for (int k = 1; k < N; k++) {
    double x = 2 * Q * R * sin(M_PI * k / N);
    // spherical Besell function of order zero 'j0' is sin(x)/x
    sph[k] = sin(x) / x;
  }
 
  std::vector<double> ratel(N);
  for (int l = 1; l < N; l++) {
    // notice that 0 < l/N < 1
    ratel[l] = rate * 4 * pow(sin(M_PI * l / N), 2);
  }
 
  const auto shift = this->getParameter("Shift");
  for (size_t i = 0; i < nData; i++) {
    double w = xValues[i] - shift;
    double S = 0.0;
    for (int l = 1; l < N; l++) {
      double lorentzian = ratel[l] / (ratel[l] * ratel[l] + w * w);
      double al = 0.0;
      for (int k = 1; k < N; k++) { // case k==N after the loop
        double y = 2 * M_PI * l * k / N;
        al += cos(y) * sph[k];
      }
      al += 1; // limit for j0 when k==N, or x==0
      al /= N;
      S += al * lorentzian;
    }
    out[i] = I * S / M_PI;
  }
}
 
void InelasticDiffRotDiscreteCircle::setWorkspace(std::shared_ptr<const API::Workspace> ws) {
  m_qValueCache.clear();
 
  auto workspace = std::dynamic_pointer_cast<const API::MatrixWorkspace>(ws);
  if (!workspace)
    return;
 
  const auto &spectrumInfo = workspace->spectrumInfo();
  const auto &detectorIDs = workspace->detectorInfo().detectorIDs();
  for (size_t idx = 0; idx < spectrumInfo.size(); idx++) {
    if (!spectrumInfo.hasDetectors(idx)) {
      m_qValueCache.clear();
      g_log.information("Cannot populate Q values from workspace - no detectors set.");
      break;
    }
 
    const auto detectorIndex = spectrumInfo.spectrumDefinition(idx)[0].first;
 
    try {
      double efixed = workspace->getEFixed(detectorIDs[detectorIndex]);
      double usingTheta = 0.5 * spectrumInfo.twoTheta(idx);
 
      double q = Mantid::Kernel::UnitConversion::convertToElasticQ(usingTheta, efixed);
 
      m_qValueCache.emplace_back(q);
    } catch (std::runtime_error &) {
      m_qValueCache.clear();
      g_log.information("Cannot populate Q values from workspace - could not "
                        "find EFixed value.");
      return;
    }
  }
}
 
} // namespace Mantid::CurveFitting::Functions