GaussianSumUtilities.h

Go to the documentation of this file.

#ifndef GaussianSumUtilities_h_
#define GaussianSumUtilities_h_

#include "TrackingTools/GsfTools/interface/MultiGaussianState1D.h"
#include "TrackingTools/GsfTools/interface/SingleGaussianState.h"
#include "TrackingTools/GsfTools/interface/MultiGaussianState.h"
#include <vector>

template <unsigned int N>
class GaussianSumUtilities {
public:
  typedef SingleGaussianState<N> SingleState;
  typedef MultiGaussianState<N> MultiState; 
//   typedef ROOT::Math::SVector<double, N> Vector;
  typedef ROOT::Math::SMatrix<double,N,N,ROOT::Math::MatRepStd<double,N> > GenMatrix;

  typedef typename SingleState::Vector Vector;
  typedef typename SingleState::Matrix Matrix;
  typedef typename MultiState::SingleStatePtr SingleStatePtr;
  typedef typename MultiState::SingleStateContainer SingleStateContainer;

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

public:
  GaussianSumUtilities (const MultiState& state) :
    theState(state), 
    theModeStatus(NotComputed) {
  } 
  ~GaussianSumUtilities () {
  }

  inline unsigned int size () const {
    return components().size();
  }
  const SingleStateContainer& components () const {
    return theState.components();
  }
  const MultiState& state () const {
    return theState;
  }
  inline double weight (unsigned int i) const {
    return components()[i]->weight();
  }
  inline const Vector& mean (unsigned int i) const {
    return components()[i]->mean();
  }
  inline const Matrix& covariance (unsigned int i) const {
    return components()[i]->covariance();
  }
  bool modeIsValid () const;
  const SingleGaussianState<N>& mode () const;
  double pdf (const Vector&) const;
  Vector d1Pdf (const Vector&) const;
  Matrix d2Pdf (const Vector&) const;
  double lnPdf (const Vector&) const;
  Vector d1LnPdf (const Vector&) const;
  Matrix d2LnPdf (const Vector&) const;

  double weight () const {
    return theState.weight();
  }
  const Vector& mean () const {
    return theState.mean();
  }
  const Matrix& covariance () const {
    return theState.covariance();
  }


protected:
  Vector computeModeWithoutTransform () const;

private:
  Matrix tensorProduct (const Vector&) const;
  double gauss (const double&, const double&, const double&) const;
  double gauss (const Vector&, 
        const Vector&,
        const Matrix&) const;
  bool findMode (Vector& mode, double& pdfAtMode,
         const Vector& xStart) const;
  void computeMode () const;
  MultiGaussianState1D constrainedState (const Vector& d,
                     const Vector& x0) const;
//   /// replacement of CLHEP determinant (which rounds off small values)
//   double determinant (const Matrix& matrix) const;
  Matrix localCovariance (const Vector& x) const;
  void setMode (const Vector& mode) const;
  void setInvalidMode () const;

  std::vector<double> pdfComponents (const Vector&) const;
  double pdf (const Vector&, const std::vector<double>&) const;
  Vector d1Pdf (const Vector&, const std::vector<double>&) const;
  Matrix d2Pdf (const Vector&, const std::vector<double>&) const;
  double lnPdf (const Vector&, const std::vector<double>&) const;
  Vector d1LnPdf (const Vector&, const std::vector<double>&) const;
  Matrix d2LnPdf (const Vector&, const std::vector<double>&) const;



private:
  const MultiState& theState;
//   int theDimension;

  mutable ModeStatus theModeStatus;
//   mutable Vector theMode;
  mutable SingleGaussianState<N> theMode;

};

#include "TrackingTools/GsfTools/interface/GaussianSumUtilities.icc"

#endif