FitPowderDiffPeaks.h

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2012 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 +
#pragma once
 
#include "MantidAPI/Algorithm.h"
#include "MantidAPI/CompositeFunction.h"
#include "MantidAPI/ITableWorkspace_fwd.h"
#include "MantidAPI/MatrixWorkspace_fwd.h"
#include "MantidCurveFitting/DllConfig.h"
#include "MantidCurveFitting/Functions/BackToBackExponential.h"
#include "MantidCurveFitting/Functions/BackgroundFunction.h"
#include "MantidCurveFitting/Functions/Polynomial.h"
#include "MantidDataObjects/TableWorkspace.h"
#include "MantidDataObjects/Workspace2D.h"
#include "MantidGeometry/Crystal/UnitCell.h"
 
#include <span>
 
namespace Mantid {
namespace CurveFitting {
namespace Algorithms {
 
class MANTID_CURVEFITTING_DLL FitPowderDiffPeaks final : public API::Algorithm {
public:
  FitPowderDiffPeaks();
 
  const std::string name() const override { return "FitPowderDiffPeaks"; }
  const std::string summary() const override { return "Fit peaks in powder diffraction pattern. "; }
 
  int version() const override { return 1; }
  const std::vector<std::string> seeAlso() const override { return {"LeBailFit"}; }
 
  const std::string category() const override { return "Diffraction\\Fitting"; }
 
private:
  void init() override;
  void exec() override;
 
  void processInputProperties();
 
  void genPeaksFromTable(const DataObjects::TableWorkspace_sptr &peakparamws);
 
  Functions::BackToBackExponential_sptr genPeak(std::map<std::string, int> hklmap,
                                                std::map<std::string, double> parammap,
                                                std::map<std::string, std::string> bk2bk2braggmap, bool &good,
                                                std::vector<int> &hkl, double &d_h);
 
  bool getHKLFromMap(std::map<std::string, int> intmap, std::vector<int> &hkl);
 
  void importInstrumentParameterFromTable(const DataObjects::TableWorkspace_sptr &parameterWS);
 
  void parseBraggPeakTable(const DataObjects::TableWorkspace_sptr &peakws,
                           std::vector<std::map<std::string, double>> &parammaps,
                           std::vector<std::map<std::string, int>> &hklmaps);
 
  void fitPeaksWithGoodStartingValues();
 
  void fitPeaksRobust();
 
  bool fitPeak(Functions::BackToBackExponential_sptr peak, Functions::BackgroundFunction_sptr background,
               double leftdev, double rightdev, size_t m_wsIndex, double &chi2);
 
  //---------------------------------------------------------------------------
 
  bool fitSinglePeakRobust(const Functions::BackToBackExponential_sptr &peak,
                           const Functions::BackgroundFunction_sptr &backgroundfunction, double peakleftbound,
                           double peakrightbound, const std::map<std::string, double> &rightpeakparammap,
                           double &finalchi2);
 
  bool fitSinglePeakSimulatedAnnealing(const Functions::BackToBackExponential_sptr &peak,
                                       const std::vector<std::string> &paramtodomc);
 
  bool fitSinglePeakConfidentX(Functions::BackToBackExponential_sptr peak);
 
  bool fitSinglePeakConfident(const Functions::BackToBackExponential_sptr &peak,
                              const Functions::BackgroundFunction_sptr &backgroundfunction, double leftbound,
                              double rightbound, double &chi2, bool &annhilatedpeak);
 
  bool fitSinglePeakConfidentY(DataObjects::Workspace2D_sptr dataws, Functions::BackToBackExponential_sptr peak,
                               double dampingfactor);
 
  bool fitOverlappedPeaks(std::vector<Functions::BackToBackExponential_sptr> peaks,
                          const Functions::BackgroundFunction_sptr &backgroundfunction, double gfwhm);
 
  bool doFitMultiplePeaks(const DataObjects::Workspace2D_sptr &dataws, size_t wsindex,
                          const API::CompositeFunction_sptr &peaksfunc,
                          std::vector<Functions::BackToBackExponential_sptr> peakfuncs, std::vector<bool> &vecfitgood,
                          std::vector<double> &vecchi2s);
 
  void estimatePeakHeightsLeBail(const DataObjects::Workspace2D_sptr &dataws, size_t wsindex,
                                 std::vector<Functions::BackToBackExponential_sptr> peaks);
 
  void setOverlappedPeaksConstraints(const std::vector<Functions::BackToBackExponential_sptr> &peaks);
 
  bool doFit1PeakSimple(const DataObjects::Workspace2D_sptr &dataws, size_t workspaceindex,
                        const Functions::BackToBackExponential_sptr &peakfunction, const std::string &minimzername,
                        size_t maxiteration, double &chi2);
 
  // bool doFit1PeakBackgroundSimple(DataObjects::Workspace2D_sptr dataws,
  // size_t workspaceindex,
  // BackToBackExponential_sptr peakfunction,
  // BackgroundFunction_sptr backgroundfunction,
  // string minimzername, size_t maxiteration, double &chi2);
 
  bool doFit1PeakBackground(const DataObjects::Workspace2D_sptr &dataws, size_t wsindex,
                            const Functions::BackToBackExponential_sptr &peak,
                            const Functions::BackgroundFunction_sptr &backgroundfunction, double &chi2);
 
  bool doFit1PeakSequential(const DataObjects::Workspace2D_sptr &dataws, size_t workspaceindex,
                            const Functions::BackToBackExponential_sptr &peakfunction,
                            std::vector<std::string> minimzernames, std::vector<size_t> maxiterations,
                            const std::vector<double> &dampfactors, double &chi2);
 
  bool doFitNPeaksSimple(const DataObjects::Workspace2D_sptr &dataws, size_t wsindex,
                         const API::CompositeFunction_sptr &peaksfunc,
                         const std::vector<Functions::BackToBackExponential_sptr> &peakfuncs,
                         const std::string &minimizername, size_t maxiteration, double &chi2);
 
  void storeFunctionParameters(const API::IFunction_sptr &function, std::map<std::string, double> &parammaps);
 
  void restoreFunctionParameters(const API::IFunction_sptr &function, std::map<std::string, double> parammap);
 
  void calculatePeakFitBoundary(size_t ileftpeak, size_t irightpeak, double &peakleftboundary,
                                double &peakrightboundary);
 
  //---------------------------------------------------------------------------
 
  bool findMaxHeight(API::MatrixWorkspace_sptr dataws, size_t wsindex, double xmin, double xmax, double &center,
                     double &centerleftbound, double &centerrightbound, int &errordirection);
 
  DataObjects::Workspace2D_sptr genPeakParameterDataWorkspace();
 
  std::pair<DataObjects::TableWorkspace_sptr, DataObjects::TableWorkspace_sptr> genPeakParametersWorkspace();
 
  void cropWorkspace(double tofmin, double tofmax);
 
  std::string parseFitParameterWorkspace(const API::ITableWorkspace_sptr &paramws);
 
  // DataObjects::Workspace2D_sptr
  // buildPartialWorkspace(API::MatrixWorkspace_sptr sourcews, size_t m_wsIndex,
  // double leftbound, double rightbound);
 
  bool estimateSinglePeakRange(Functions::BackToBackExponential_sptr peak,
                               Functions::BackgroundFunction_sptr background,
                               Functions::BackToBackExponential_sptr rightpeak, double fwhm, bool ismostright,
                               size_t m_wsIndex, double &chi2);
 
  void observePeakRange(const Functions::BackToBackExponential_sptr &thispeak,
                        const Functions::BackToBackExponential_sptr &rightpeak, double refpeakshift,
                        double &peakleftbound, double &peakrightbound);
 
  // void estimateBackground(DataObjects::Workspace2D_sptr dataws);
 
  void subtractBackground(DataObjects::Workspace2D_sptr dataws);
 
  bool estimateFWHM(DataObjects::Workspace2D_sptr dataws, size_t wsindex, double tof_h, double &leftfwhm,
                    double &rightfwhm);
 
  bool doFitBackground(DataObjects::Workspace2D_sptr dataws, Functions::BackgroundFunction_sptr background,
                       double leftpeakbound, double rightpeakbound);
 
  std::pair<bool, double> doFitPeak_Old(DataObjects::Workspace2D_sptr dataws,
                                        Functions::BackToBackExponential_sptr peak, double guessedfwhm, bool calchi2);
 
  std::pair<bool, double> doFitPeak(const DataObjects::Workspace2D_sptr &dataws,
                                    const Functions::BackToBackExponential_sptr &peakfunction, double guessedfwhm);
 
  bool doFitGaussianPeak(const DataObjects::Workspace2D_sptr &dataws, size_t workspaceindex, double in_center,
                         double leftfwhm, double rightfwhm, double &center, double &sigma, double &height);
 
  DataObjects::Workspace2D_sptr genOutputFittedPatternWorkspace(const std::vector<double> &pattern, int workspaceindex);
 
  void calculate1PeakGroup(std::vector<size_t> peakindexes, Functions::BackgroundFunction_sptr background);
 
  std::string parseFitResult(const API::IAlgorithm_sptr &fitalg, double &chi2, bool &fitsuccess);
 
  double calculatePeakCentreTOF(int h, int k, int l);
 
  double getParameter(const std::string &parname);
 
  void fitPeaksGroup(std::vector<size_t> peakindexes);
 
  DataObjects::Workspace2D_sptr buildPartialWorkspace(const API::MatrixWorkspace_sptr &sourcews, size_t workspaceindex,
                                                      double leftbound, double rightbound);
 
  void plotFunction(const API::IFunction_sptr &peakfunction, const Functions::BackgroundFunction_sptr &background,
                    const API::FunctionDomain1DVector &domain);
 
  //-----------------------------------------------------------------------------------------------
 
  API::MatrixWorkspace_sptr m_dataWS;
 
  DataObjects::TableWorkspace_sptr m_peakParamTable;
 
  DataObjects::TableWorkspace_sptr m_profileTable;
 
  // Map for all peaks to fit individually
  // Disabled std::map<std::vector<int>,
  // CurveFitting::BackToBackExponential_sptr> m_peaksmap;
 
  std::vector<std::pair<double, std::pair<std::vector<int>, Functions::BackToBackExponential_sptr>>> m_vecPeakFunctions;
 
  std::vector<double> m_peakFitChi2;
 
  std::vector<bool> m_goodFit;
 
  std::map<std::string, double> m_instrumentParmaeters;
 
  std::vector<double> m_peakData;
 
  std::vector<std::string> mPeakParameterNames;
 
  int m_wsIndex;
 
  double m_tofMin;
  double m_tofMax;
 
  bool m_useGivenTOFh;
 
  bool m_confidentInInstrumentParameters;
 
  std::vector<int> m_minimumHKL;
 
  int m_numPeaksLowerToMin;
 
  std::vector<size_t> m_indexGoodFitPeaks;
  std::vector<double> m_chi2GoodFitPeaks;
 
  enum { ROBUSTFIT, TRUSTINPUTFIT } m_fitMode;
 
  enum { HKLCALCULATION, FROMBRAGGTABLE } m_genPeakStartingValue;
 
  std::vector<int> m_rightmostPeakHKL;
 
  double m_rightmostPeakLeftBound;
 
  double m_rightmostPeakRightBound;
 
  double m_minPeakHeight;
 
  Geometry::UnitCell m_unitCell;
 
  bool m_fitPeakBackgroundComposite;
};
 
inline double linearInterpolateX(double x0, double xf, double y0, double yf, double y) {
  double x = ((xf - x0) * y - (xf * y0 - x0 * yf)) / (yf - y0);
  return x;
}
 
inline double linearInterpolateY(double x0, double xf, double y0, double yf, double x) {
  double y = ((xf * y0 - x0 * yf) + x * (yf - y0)) / (xf - x0);
  return y;
}
 
void estimateBackgroundCoarse(const DataObjects::Workspace2D_sptr &dataws,
                              const Functions::BackgroundFunction_sptr &background, size_t wsindexraw,
                              size_t wsindexbkgd, size_t wsindexpeak);
 
bool observePeakParameters(const DataObjects::Workspace2D_sptr &dataws, size_t wsindex, double &centre, double &height,
                           double &fwhm, std::string &errmsg);
 
size_t findMaxValue(std::span<double const> Y);
 
size_t findMaxValue(const API::MatrixWorkspace_sptr &dataws, size_t wsindex, double leftbound, double rightbound);
 
std::string getFunctionInfo(const API::IFunction_sptr &function);
 
} // namespace Algorithms
} // namespace CurveFitting
} // namespace Mantid