FitPeaks.h

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 +
#pragma once
 
#include "MantidAPI/Algorithm.h"
#include "MantidAPI/IBackgroundFunction.h"
#include "MantidAPI/IPeakFunction.h"
#include "MantidAPI/ITableWorkspace.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/MatrixWorkspace_fwd.h"
#include "MantidAPI/MultiDomainFunction.h"
#include "MantidAlgorithms/DllConfig.h"
#include "MantidAlgorithms/PeakParameterHelper.h"
#include "MantidDataObjects/EventWorkspace.h"
#include "MantidDataObjects/TableWorkspace.h"
#include "MantidKernel/cow_ptr.h"
 
#include <utility>
 
namespace Mantid {
namespace HistogramData {
class HistogramX;
class HistogramY;
} // namespace HistogramData
 
namespace Algorithms {
 
namespace FitPeaksAlgorithm {
struct FitFunction {
  API::IPeakFunction_sptr peakfunction;
  API::IBackgroundFunction_sptr bkgdfunction;
};
 
class PeakFitResult {
public:
  PeakFitResult(size_t num_peaks, size_t num_params);
  double getPeakPosition(size_t ipeak) const;
  double getCost(size_t ipeak) const;
  size_t getNumberParameters() const;
  size_t getNumberPeaks() const;
  double getParameterValue(size_t ipeak, size_t iparam) const;
  double getParameterError(size_t ipeak, size_t iparam) const;
  void setRecord(size_t ipeak, const double cost, const double peak_position, const FitFunction &fit_functions);
  void setBadRecord(size_t ipeak, const double peak_position);
  void setFunctionParameters(size_t ipeak, std::vector<double> &param_values);
 
private:
  size_t m_function_parameters_number;
  // goodness of fitting
  std::vector<double> m_costs;
  // fitted peak positions
  std::vector<double> m_fitted_peak_positions;
  // fitted peak and background parameters
  std::vector<std::vector<double>> m_function_parameters_vector;
  std::vector<std::vector<double>> m_function_errors_vector;
};
 
class PeakFitPreCheckResult {
public:
  PeakFitPreCheckResult()
      : m_submitted_spectrum_peaks{0}, m_submitted_individual_peaks{0}, m_low_count_spectrum{0}, m_out_of_range{0},
        m_low_count_individual{0}, m_not_enough_datapoints{0}, m_low_snr{0} {}
  PeakFitPreCheckResult &operator+=(const PeakFitPreCheckResult &another);
 
public:
  void setNumberOfSubmittedSpectrumPeaks(const size_t n);
  void setNumberOfSubmittedIndividualPeaks(const size_t n);
  void setNumberOfSpectrumPeaksWithLowCount(const size_t n);
  void setNumberOfOutOfRangePeaks(const size_t n);
  void setNumberOfIndividualPeaksWithLowCount(const size_t n);
  void setNumberOfPeaksWithNotEnoughDataPoints(const size_t n);
  void setNumberOfPeaksWithLowSignalToNoise(const size_t n);
  bool isIndividualPeakRejected() const;
  std::string getReport() const;
 
private:
  // number of peaks submitted for spectrum fitting
  size_t m_submitted_spectrum_peaks;
  // number of peaks submitted for individual fitting. Since some spectra might fail a pre-check, not all peaks might
  // make it to the individual fitting
  size_t m_submitted_individual_peaks;
  // number of peaks rejected as a whole spectrum due to its low signal count
  size_t m_low_count_spectrum;
  // number of peaks rejected individually because their predicted position is out of range
  size_t m_out_of_range;
  // number of peaks rejected individually due to low signal count
  size_t m_low_count_individual;
  // number of peask rejected due to not enough data points
  size_t m_not_enough_datapoints;
  // number of peaks rejected due to low signal-to-noise ratio
  size_t m_low_snr;
};
} // namespace FitPeaksAlgorithm
 
class MANTID_ALGORITHMS_DLL FitPeaks final : public API::Algorithm {
public:
  FitPeaks();
 
  const std::string name() const override { return "FitPeaks"; }
 
  const std::string summary() const override { return "Fit one or multiple peaks in all spectra of a given workspace"; }
 
  int version() const override { return (1); }
 
  const std::string category() const override { return "Optimization"; }
 
  std::map<std::string, std::string> validateInputs() override;
 
  static bool fitStatusIsConverged(const std::string &fitStatus, const bool strict);
 
private:
  void init() override;
  void exec() override;
 
  void processInputs();
  void processInputPeakCenters();
  void processInputPeakTolerance();
  void processInputFunctions();
  void processInputFitRanges();
 
  void generateFittedParametersValueWorkspaces();
  void generateOutputPeakPositionWS();
  void generateCalculatedPeaksWS();
 
  void convertParametersNameToIndex();
 
  std::vector<std::shared_ptr<FitPeaksAlgorithm::PeakFitResult>> fitPeaks();
 
  void fitSpectrumPeaks(size_t wi, const std::vector<double> &expected_peak_centers,
                        const std::shared_ptr<FitPeaksAlgorithm::PeakFitResult> &fit_result,
                        std::vector<std::vector<double>> &lastGoodPeakParameters,
                        std::vector<size_t> &lastGoodPeakSpectra,
                        const std::shared_ptr<FitPeaksAlgorithm::PeakFitPreCheckResult> &pre_check_result);
 
  bool fitBackground(const size_t &ws_index, const std::pair<double, double> &fit_window,
                     const double &expected_peak_pos, const API::IBackgroundFunction_sptr &bkgd_func);
 
  // Peak fitting suite
  double fitIndividualPeak(size_t wi, const API::IAlgorithm_sptr &fitter, const double expected_peak_center,
                           const double peak_pos_tolerance, const std::pair<double, double> &fitwindow,
                           const bool estimate_peak_width, const API::IPeakFunction_sptr &peakfunction,
                           const API::IBackgroundFunction_sptr &bkgdfunc,
                           const std::shared_ptr<FitPeaksAlgorithm::PeakFitPreCheckResult> &pre_check_result);
 
  double fitFunctionSD(const API::IAlgorithm_sptr &fit, const API::IPeakFunction_sptr &peak_function,
                       const API::IBackgroundFunction_sptr &bkgd_function, const API::MatrixWorkspace_sptr &dataws,
                       size_t wsindex, const std::pair<double, double> &peak_range, const double &expected_peak_center,
                       const double peak_pos_tolerance, bool estimate_peak_width, bool estimate_background);
 
  void recalculateErrorsWithoutConstraint(const API::IPeakFunction_sptr &peak_function,
                                          const API::IBackgroundFunction_sptr &bkgd_function,
                                          const API::MatrixWorkspace_sptr &dataws, size_t wsindex,
                                          const std::pair<double, double> &peak_range);
 
  double fitFunctionMD(API::IFunction_sptr fit_function, const API::MatrixWorkspace_sptr &dataws, const size_t wsindex,
                       const std::pair<double, double> &vec_xmin, const std::pair<double, double> &vec_xmax);
 
  double fitFunctionHighBackground(const API::IAlgorithm_sptr &fit, const std::pair<double, double> &fit_window,
                                   const size_t &ws_index, const double &expected_peak_center,
                                   const double peak_pos_tolerance, bool observe_peak_shape,
                                   const API::IPeakFunction_sptr &peakfunction,
                                   const API::IBackgroundFunction_sptr &bkgdfunc);
 
  void setupParameterTableWorkspace(const API::ITableWorkspace_sptr &table_ws,
                                    const std::vector<std::string> &param_names, bool with_chi2);
 
  void histRangeToIndexBounds(size_t iws, const std::pair<double, double> &range, size_t &left_index,
                              size_t &right_index); 
 
  size_t histRangeToDataPointCount(size_t iws, const std::pair<double, double> &range);
 
  void getRangeData(size_t iws, const std::pair<double, double> &range, std::vector<double> &vec_x,
                    std::vector<double> &vec_y, std::vector<double> &vec_e);
 
  double numberCounts(size_t iws);
 
  double numberCounts(size_t iws, const std::pair<double, double> &range);
 
  double calculateSignalToNoiseRatio(size_t iws, const std::pair<double, double> &range,
                                     const API::IBackgroundFunction_sptr &bkgd_function);
 
  API::MatrixWorkspace_sptr createMatrixWorkspace(const std::vector<double> &vec_x, const std::vector<double> &vec_y,
                                                  const std::vector<double> &vec_e);
 
  bool decideToEstimatePeakParams(const bool firstPeakInSpectrum, const size_t wsindex,
                                  const API::IPeakFunction_sptr &peak_function);
 
  bool processSinglePeakFitResult(size_t wsindex, size_t peakindex, const double cost,
                                  const std::vector<double> &expected_peak_positions,
                                  const FitPeaksAlgorithm::FitFunction &fitfunction,
                                  const std::shared_ptr<FitPeaksAlgorithm::PeakFitResult> &fit_result);
 
  void calculateFittedPeaks(const std::vector<std::shared_ptr<FitPeaksAlgorithm::PeakFitResult>> &fit_results);
 
  double calculateSignalToSigmaRatio(const size_t &iws, const std::pair<double, double> &peakWindow,
                                     const API::IPeakFunction_sptr &peakFunction);
 
  std::string getPeakHeightParameterName(const API::IPeakFunction_const_sptr &peak_function);
 
  void processOutputs(std::vector<std::shared_ptr<FitPeaksAlgorithm::PeakFitResult>> fit_result_vec);
 
  void writeFitResult(size_t wi, const std::vector<double> &expected_positions,
                      const std::shared_ptr<FitPeaksAlgorithm::PeakFitResult> &fit_result);
 
  bool isObservablePeakProfile(const std::string &peakprofile);
 
  // log a message disregarding the current logging offset
  void logNoOffset(const size_t &priority, const std::string &msg);
 
  // create a Fit child alg
  API::IAlgorithm_sptr createChildFit();
 
  //------- Workspaces-------------------------------------
  API::MatrixWorkspace_sptr m_inputMatrixWS;
  bool m_inputIsDSpace;
  DataObjects::EventWorkspace_const_sptr m_inputEventWS; // cast from m_inputWS
  API::MatrixWorkspace_sptr m_outputPeakPositionWorkspace; // output workspace for peak positions
  API::ITableWorkspace_sptr m_fittedParamTable;
  API::ITableWorkspace_sptr m_fitErrorTable;
  bool m_rawPeaksTable;
  API::MatrixWorkspace_sptr m_fittedPeakWS;
 
  //-------- Functions ------------------------------------------------------
  API::IPeakFunction_sptr m_peakFunction;
  API::IBackgroundFunction_sptr m_bkgdFunction;
  API::IBackgroundFunction_sptr m_linearBackgroundFunction;
 
  std::string m_minimizer;
  std::string m_costFunction;
  bool m_strictConvergence{true};
  bool m_fitPeaksFromRight;
  int m_fitIterations;
  // Copy the peak parameters from the last successfully fit peak
  bool m_copyLastGoodPeakParameters;
  // If true, parameters marked as fixed on the peak function (e.g. parameters
  // calculated from the instrument geometry) remain fixed during fitting.
  // If false (default), such parameters are unfixed so they can be refined.
  bool m_respectFixedPeakParameters;
 
  //-------- Input param init values --------------------------------
  std::vector<size_t> m_initParamIndexes;
 
  std::vector<double> m_peakCenters;
  API::MatrixWorkspace_const_sptr m_peakCenterWorkspace;
  std::size_t m_numPeaksToFit;
  bool m_uniformPeakPositions;
 
  std::function<std::vector<double>(std::size_t const &)> m_getExpectedPeakPositions;
  std::function<std::pair<double, double>(std::size_t const &, std::size_t const &)> m_getPeakFitWindow;
  void checkWorkspaceIndices(std::size_t const &);
  void checkPeakIndices(std::size_t const &, std::size_t const &);
  void checkPeakWindowEdgeOrder(double const &, double const &);
 
  double m_peakWidthPercentage;
 
  //--------- Fitting range -----------------------------------------
  std::size_t m_startWorkspaceIndex;
  std::size_t m_stopWorkspaceIndex;
  std::size_t m_numSpectraToFit;
  std::vector<double> m_peakPosTolerances;
 
  Algorithms::PeakParameterHelper::EstimatePeakWidth m_peakWidthEstimateApproach;
  bool m_constrainPeaksPosition;
  bool m_constrainByPositionTolerance{false};
  bool m_fractionalPositionTolerance{false};
  bool m_calculateUnconstrainedErrors{false};
 
  std::vector<std::vector<double>> m_peakWindowVector;
  API::MatrixWorkspace_const_sptr m_peakWindowWorkspace;
 
  std::vector<std::string> m_peakParamNames;
  std::vector<double> m_initParamValues;
  API::ITableWorkspace_const_sptr m_profileStartingValueTable;
  bool m_uniformProfileStartingValue;
 
  // Criteria for fitting peaks
  double m_minPeakHeight;
 
  // Criteria for rejecting non-peaks or weak peaks from fitting
  double m_minSignalToNoiseRatio;
  double m_minPeakTotalCount;
 
  double m_minSignalToSigmaRatio;
 
  bool m_highBackground;
 
  //----- Result criterias ---------------
  bool m_peakPosTolCase234;
};
 
} // namespace Algorithms
} // namespace Mantid