CrystalFieldPeakUtils.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/CrystalFieldPeakUtils.h"
 
#include "MantidAPI/CompositeFunction.h"
#include "MantidAPI/FunctionFactory.h"
#include "MantidAPI/IPeakFunction.h"
#include "MantidCurveFitting/Constraints/BoundaryConstraint.h"
 
#include <algorithm>
#include <cmath>
 
namespace Mantid::CurveFitting::Functions::CrystalFieldUtils {
 
using namespace CurveFitting;
using namespace Kernel;
using namespace API;
 
double calculateWidth(double x, const std::vector<double> &xVec, const std::vector<double> &yVec) {
  assert(xVec.size() == yVec.size());
  auto upperIt = std::lower_bound(xVec.begin(), xVec.end(), x);
  if (upperIt == xVec.end() || x < xVec.front()) {
    return -1.0;
  }
  if (upperIt == xVec.begin()) {
    return yVec.front();
  }
  double lowerX = *(upperIt - 1);
  double upperX = *upperIt;
  auto i = std::distance(xVec.begin(), upperIt) - 1;
  return yVec[i] + (yVec[i + 1] - yVec[i]) / (upperX - lowerX) * (x - lowerX);
}
 
void setWidthConstraint(API::IPeakFunction &peak, double fwhm, double fwhmVariation) {
  double upperBound = fwhm + fwhmVariation;
  double lowerBound = fwhm - fwhmVariation;
  bool fix = lowerBound == upperBound;
  if (!fix) {
    if (lowerBound < 0.0) {
      lowerBound = 0.0;
    }
    if (lowerBound >= upperBound) {
      lowerBound = upperBound / 2;
    }
  }
  if (peak.name() == "Lorentzian") {
    if (fix) {
      peak.fixParameter("FWHM");
      return;
    }
    peak.removeConstraint("FWHM");
    auto constraint = std::make_unique<Constraints::BoundaryConstraint>(&peak, "FWHM", lowerBound, upperBound);
    peak.addConstraint(std::move(constraint));
  } else if (peak.name() == "Gaussian") {
    if (fix) {
      peak.fixParameter("Sigma");
      return;
    }
    const double WIDTH_TO_SIGMA = 2.0 * sqrt(2.0 * M_LN2);
    lowerBound /= WIDTH_TO_SIGMA;
    upperBound /= WIDTH_TO_SIGMA;
    peak.removeConstraint("Sigma");
    auto constraint = std::make_unique<Constraints::BoundaryConstraint>(&peak, "Sigma", lowerBound, upperBound);
    peak.addConstraint(std::move(constraint));
  } else {
    throw std::runtime_error("Cannot set constraint on width of " + peak.name());
  }
}
 
size_t calculateNPeaks(const API::FunctionValues &centresAndIntensities) { return centresAndIntensities.size() / 2; }
 
size_t calculateMaxNPeaks(size_t nPeaks) { return nPeaks + nPeaks / 2 + 1; }
 
inline void ignorePeak(API::IPeakFunction &peak, double fwhm) {
  peak.setHeight(0.0);
  peak.fixAll(true);
  peak.setFwhm(fwhm);
}
 
void setPeakWidth(API::IPeakFunction &peak, double centre, const std::vector<double> &xVec,
                  const std::vector<double> &yVec, double fwhmVariation, double defaultFWHM, bool useDefaultFWHM) {
  if (useDefaultFWHM) {
    peak.setFwhm(defaultFWHM);
  } else {
    auto fwhm = calculateWidth(centre, xVec, yVec);
    if (fwhm > 0.0) {
      peak.setFwhm(fwhm);
      setWidthConstraint(peak, fwhm, fwhmVariation);
    } else {
      ignorePeak(peak, defaultFWHM);
    }
  }
}
 
API::IPeakFunction_sptr createPeak(const std::string &peakShape, double centre, double intensity,
                                   const std::vector<double> &xVec, const std::vector<double> &yVec,
                                   double fwhmVariation, double defaultFWHM, bool isGood, bool fixAllPeaks) {
  auto fun = API::FunctionFactory::Instance().createFunction(peakShape);
  auto peak = std::dynamic_pointer_cast<API::IPeakFunction>(fun);
  if (!peak) {
    throw std::runtime_error("A peak function is expected.");
  }
  bool useDefaultFWHM = xVec.empty();
  const bool fixByDefault = true;
  if (isGood) {
    peak->setCentre(centre);
    peak->setIntensity(intensity);
    setPeakWidth(*peak, centre, xVec, yVec, fwhmVariation, defaultFWHM, useDefaultFWHM);
    peak->fixCentre(fixByDefault);
    peak->fixIntensity(fixByDefault);
  } else {
    ignorePeak(*peak, defaultFWHM);
  }
  if (fixAllPeaks) {
    peak->fixAll(fixByDefault);
  }
  return peak;
}
 
size_t buildSpectrumFunction(API::CompositeFunction &spectrum, const std::string &peakShape,
                             const API::FunctionValues &centresAndIntensities, const std::vector<double> &xVec,
                             const std::vector<double> &yVec, double fwhmVariation, double defaultFWHM,
                             size_t nRequiredPeaks, bool fixAllPeaks) {
  if (xVec.size() != yVec.size()) {
    throw std::runtime_error("WidthX and WidthY must have the same size.");
  }
 
  auto nPeaks = calculateNPeaks(centresAndIntensities);
  auto maxNPeaks = calculateMaxNPeaks(nPeaks);
  if (nRequiredPeaks > maxNPeaks) {
    maxNPeaks = nRequiredPeaks;
  }
  for (size_t i = 0; i < maxNPeaks; ++i) {
    const bool isGood = i < nPeaks;
    const auto centre = isGood ? centresAndIntensities.getCalculated(i) : 0.0;
    const auto intensity = isGood ? centresAndIntensities.getCalculated(i + nPeaks) : 0.0;
    auto peak = createPeak(peakShape, centre, intensity, xVec, yVec, fwhmVariation, defaultFWHM, isGood, fixAllPeaks);
    spectrum.addFunction(peak);
  }
  return nPeaks;
}
 
void updatePeakWidth(API::IPeakFunction &peak, double centre, const std::vector<double> &xVec,
                     const std::vector<double> &yVec, double fwhmVariation) {
  bool mustUpdateWidth = !xVec.empty();
  if (mustUpdateWidth) {
    auto fwhm = peak.fwhm();
    auto expectedFwhm = calculateWidth(centre, xVec, yVec);
    if (expectedFwhm <= 0.0) {
      ignorePeak(peak, fwhm);
    } else if (fabs(fwhm - expectedFwhm) > fwhmVariation) {
      peak.setFwhm(expectedFwhm);
      setWidthConstraint(peak, expectedFwhm, fwhmVariation);
    }
  }
}
 
void updatePeak(API::IPeakFunction &peak, double centre, double intensity, const std::vector<double> &xVec,
                const std::vector<double> &yVec, double fwhmVariation, bool isGood, bool fixAllPeaks) {
  const bool fixByDefault = true;
  if (isGood) {
    peak.unfixAllDefault();
    peak.setCentre(centre);
    peak.setIntensity(intensity);
    updatePeakWidth(peak, centre, xVec, yVec, fwhmVariation);
    peak.unfixIntensity();
    peak.fixIntensity(fixByDefault);
    peak.fixCentre(fixByDefault);
    if (fixAllPeaks) {
      peak.fixAll(fixByDefault);
    }
  } else {
    peak.setHeight(0.0);
    peak.fixIntensity(fixByDefault);
    peak.fixCentre(fixByDefault);
  }
}
 
size_t updateSpectrumFunction(API::CompositeFunction &spectrum, const std::string &peakShape,
                              const FunctionValues &centresAndIntensities, size_t iFirst,
                              const std::vector<double> &xVec, const std::vector<double> &yVec, double fwhmVariation,
                              double defaultFWHM, bool fixAllPeaks) {
  size_t nGoodPeaks = calculateNPeaks(centresAndIntensities);
  size_t maxNPeaks = calculateMaxNPeaks(nGoodPeaks);
  size_t nFunctions = spectrum.nFunctions();
 
  for (size_t i = 0; i < maxNPeaks; ++i) {
    const bool isGood = i < nGoodPeaks;
    auto centre = isGood ? centresAndIntensities.getCalculated(i) : 0.0;
    auto intensity = isGood ? centresAndIntensities.getCalculated(i + nGoodPeaks) : 0.0;
 
    if (i < nFunctions - iFirst) {
      auto fun = spectrum.getFunction(i + iFirst);
      auto &peak = dynamic_cast<API::IPeakFunction &>(*fun);
      updatePeak(peak, centre, intensity, xVec, yVec, fwhmVariation, isGood, fixAllPeaks);
    } else {
      auto peakPtr =
          createPeak(peakShape, centre, intensity, xVec, yVec, fwhmVariation, defaultFWHM, isGood, fixAllPeaks);
      spectrum.addFunction(peakPtr);
    }
  }
  // If there are any peaks above the maxNPeaks, ignore them
  // but don't remove
  const bool fixByDefault = true;
  for (size_t i = maxNPeaks; i < nFunctions - iFirst; ++i) {
    auto fun = spectrum.getFunction(i + iFirst);
    auto &peak = dynamic_cast<API::IPeakFunction &>(*fun);
    const auto fwhm = peak.fwhm();
    peak.setHeight(0.0);
    peak.setFwhm(fwhm);
    peak.fixIntensity(fixByDefault);
    peak.fixCentre(fixByDefault);
  }
  return nGoodPeaks;
}
 
} // namespace Mantid::CurveFitting::Functions::CrystalFieldUtils