MultiGaussianStateCombiner1D Class Reference

#include <MultiGaussianStateCombiner1D.h>

Public Member Functions
SingleGaussianState1D	combine (const MultiGaussianState1D &theState) const
SingleGaussianState1D	combine (const VSC &theComponents) const

Private Types
typedef std::vector < SingleGaussianState1D >	VSC

Detailed Description

Class to collapse (combine) a Gaussian mixture of states into one. (c.f. R. Fruewirth et.al., Comp.Phys.Comm 100 (1997) 1

Definition at line 13 of file MultiGaussianStateCombiner1D.h.

Member Typedef Documentation

typedef std::vector<SingleGaussianState1D> MultiGaussianStateCombiner1D::VSC

private

Definition at line 16 of file MultiGaussianStateCombiner1D.h.

Member Function Documentation

SingleGaussianState1D MultiGaussianStateCombiner1D::combine

(

const MultiGaussianState1D &

theState

)

const

Definition at line 9 of file MultiGaussianStateCombiner1D.cc.

References MultiGaussianState1D::components().

Referenced by MultiGaussianState1D::checkCombinedState().

{
  return combine(theState.components());
}

SingleGaussianState1D MultiGaussianStateCombiner1D::combine

(

const VSC &

theComponents

)

const

Definition at line 15 of file MultiGaussianStateCombiner1D.cc.

References gather_cfg::cout, Exception, and puppiForMET_cff::weight.

{
  if (theComponents.empty()) {
    throw cms::Exception("LogicError") 
      << "MultiGaussianStateCombiner1D:: state container to collapse is empty";
//     return SingleState(SingleState::Vector(),SingleState::Matrix(),0.);
//     return SingleState(SingleState::Vector(),
//             SingleState::Matrix(),0.);
    return SingleGaussianState1D();
  }

  const SingleGaussianState1D firstState(theComponents.front());
  if (theComponents.size()==1)  return firstState;

//   int size = firstState.mean().num_row();
  double meanMean(0.);
  double weightSum(0.);
//   double weight;
  double measCovar1(0.);
  double measCovar2(0.);
  for ( VSC::const_iterator mixtureIter1 = theComponents.begin();
    mixtureIter1 != theComponents.end(); mixtureIter1++ ) {
    double weight = mixtureIter1->weight();
    weightSum += weight;

    double mean1 = mixtureIter1->mean();
    meanMean += weight * mean1;
    measCovar1 += weight * mixtureIter1->variance();

    for ( VSC::const_iterator mixtureIter2 = mixtureIter1+1;
    mixtureIter2 != theComponents.end(); mixtureIter2++ ) {
      double posDiff = mean1 - mixtureIter2->mean();
//       SingleState::Matrix s(1,1); //stupid trick to make CLHEP work decently
//       measCovar2 +=weight * (*mixtureIter2).weight() *
//                      s.similarity(posDiff.T().T());
//       SingleState::Matrix mat;
//       for ( unsigned int i1=0; i1<N; i1++ ) {
//  for ( unsigned int i2=0; i2<=i1; i2++ )  mat(i1,i2) = posDiff(i1)*posDiff(i2);
//       }
//       measCovar2 += weight * (*mixtureIter2).weight() * mat;
      // 
      // TensorProd yields a general matrix - need to convert to a symm. matrix
      double covGen = posDiff*posDiff;
//       double covSym(covGen.LowerBlock());
      measCovar2 += weight * mixtureIter2->weight() * covGen;
    }
  }

  double measCovar;
  if (weightSum<DBL_MIN){
    std::cout << "MultiGaussianStateCombiner1D:: New state has total weight of 0."
          << std::endl;
//     meanMean = SingleState::Vector(size,0);
    meanMean = 0.;
    measCovar = 0.;
    weightSum = 0.;
  } else {
    meanMean /= weightSum;
    measCovar1 *= (1./weightSum);
    measCovar2 *= (1./weightSum/weightSum);
    measCovar = measCovar1 + measCovar2;
//     measCovar = measCovar1/weightSum + measCovar2/weightSum/weightSum;
  }

  return SingleGaussianState1D(meanMean, measCovar, weightSum);
}