---

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>

// #define DRAW_GSND

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) {
#ifdef DRAW_GSND
//     v1Draw_ = 1;
//     v2Draw_ = 2;
    instance_ = this;
#endif
  } 
  ~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;


#ifdef DRAW_GSND
public:
  void setDraw (int i1, int i2) {v1Draw_=i1; v2Draw_=i2;}
private:
  static double fcn2 (double* x, double* p);
  static double fcn2mode (double* x, double* p);
#endif

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

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

#ifdef DRAW_GSND
  static int v1Draw_;
  static int v2Draw_;
  static GaussianSumUtilities<N>* instance_;
#endif
};

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

#endif

---