GaussianSumUtilities1D.h

Go to the documentation of this file.

#ifndef GaussianSumUtilities1D_h_
#define GaussianSumUtilities1D_h_

// #include "TROOT.h"

#include "TrackingTools/GsfTools/interface/SingleGaussianState1D.h"
#include "TrackingTools/GsfTools/interface/MultiGaussianState1D.h"

#include <vector>

class GaussianSumUtilities1D {
private:
  enum ModeStatus { Valid, NotValid, NotComputed };

public:
  GaussianSumUtilities1D(const MultiGaussianState1D& state)
      : theState(state),
        //     theStates(state.components()),
        theModeStatus(NotComputed) {}
  ~GaussianSumUtilities1D() {}

  inline unsigned int size() const { return components().size(); }
  inline const std::vector<SingleGaussianState1D>& components() const { return theState.components(); }
  inline double weight(unsigned int i) const { return components()[i].weight(); }
  inline double mean(unsigned int i) const { return components()[i].mean(); }
  inline double standardDeviation(unsigned int i) const {
    //     return sqrt(components()[i].variance());
    return components()[i].standardDeviation();
  }
  inline double variance(unsigned int i) const { return components()[i].variance(); }
  double pdf(unsigned int i, double x) const;
  double quantile(const double) const;
  bool modeIsValid() const;
  const SingleGaussianState1D& mode() const;
  double pdf(double) const;
  double cdf(const double&) const;
  double d1Pdf(const double&) const;
  double d2Pdf(const double&) const;
  double d3Pdf(const double&) const;
  double lnPdf(const double&) const;
  double d1LnPdf(const double&) const;
  double d2LnPdf(const double&) const;

  double weight() const { return theState.weight(); }
  double mean() const { return theState.mean(); }
  double variance() const { return theState.variance(); }

private:
  bool findMode(double& mode, double& pdfAtMode, const double& xStart, const double& scale) const;
  static double gauss(double, double, double);
  static double gaussInt(double, double, double);
  double combinedMean() const;
  void computeMode() const;
  double localVariance(double x) const;

  // the state of the mode finder
  struct FinderState {
    FinderState() {}
    FinderState(size_t n) : pdfs(n) {}
    double x;
    double y;
    double yd;   // d1LnPdf
    double yd2;  // d2LnPdf
    std::vector<double> pdfs;
  };

  // update tre state at x
  void update(FinderState& state, double x) const;

  std::vector<double> pdfComponents(const double&) const;
  void pdfComponents(double, std::vector<double>&) const;
  static double pdf(double, const std::vector<double>&);
  double d1Pdf(double, const std::vector<double>&) const;
  double d2Pdf(double, const std::vector<double>&) const;
  double d3Pdf(double, const std::vector<double>&) const;
  static double lnPdf(double, const std::vector<double>&);
  double d1LnPdf(double, const std::vector<double>&) const;
  double d2LnPdf(double, const std::vector<double>&) const;

private:
  const MultiGaussianState1D& theState;
  //   std::vector<SingleGaussianState1D> theStates;

  mutable ModeStatus theModeStatus;
  mutable SingleGaussianState1D theMode;
  //   mutable double theMode;
};
#endif