GsfTools/src/MultiGaussianStateCombiner1D.cc

0001 #include "TrackingTools/GsfTools/interface/MultiGaussianStateCombiner1D.h"
0002
0003 #include "FWCore/Utilities/interface/Exception.h"
0004
0005 #include <cfloat>
0006
0007 SingleGaussianState1D MultiGaussianStateCombiner1D::combine(const MultiGaussianState1D& theState) const {
0008   return combine(theState.components());
0009 }
0010
0011 SingleGaussianState1D MultiGaussianStateCombiner1D::combine(const VSC& theComponents) const {
0012   if (theComponents.empty()) {
0013     throw cms::Exception("LogicError") << "MultiGaussianStateCombiner1D:: state container to collapse is empty";
0014     //     return SingleState(SingleState::Vector(),SingleState::Matrix(),0.);
0015     //     return SingleState(SingleState::Vector(),
0016     //             SingleState::Matrix(),0.);
0017     return SingleGaussianState1D();
0018   }
0019
0020   const SingleGaussianState1D& firstState(theComponents.front());
0021   if (theComponents.size() == 1)
0022     return firstState;
0023
0024   //   int size = firstState.mean().num_row();
0025   double meanMean(0.);
0026   double weightSum(0.);
0027   //   double weight;
0028   double measCovar1(0.);
0029   double measCovar2(0.);
0030   for (VSC::const_iterator mixtureIter1 = theComponents.begin(); mixtureIter1 != theComponents.end(); mixtureIter1++) {
0031     double weight = mixtureIter1->weight();
0032     weightSum += weight;
0033
0034     double mean1 = mixtureIter1->mean();
0035     meanMean += weight * mean1;
0036     measCovar1 += weight * mixtureIter1->variance();
0037
0038     for (VSC::const_iterator mixtureIter2 = mixtureIter1 + 1; mixtureIter2 != theComponents.end(); mixtureIter2++) {
0039       double posDiff = mean1 - mixtureIter2->mean();
0040       //       SingleState::Matrix s(1,1); //stupid trick to make CLHEP work decently
0041       //       measCovar2 +=weight * (*mixtureIter2).weight() *
0042       //                    s.similarity(posDiff.T().T());
0043       //       SingleState::Matrix mat;
0044       //       for ( unsigned int i1=0; i1<N; i1++ ) {
0045       //    for ( unsigned int i2=0; i2<=i1; i2++ )  mat(i1,i2) = posDiff(i1)*posDiff(i2);
0046       //       }
0047       //       measCovar2 += weight * (*mixtureIter2).weight() * mat;
0048       //
0049       // TensorProd yields a general matrix - need to convert to a symm. matrix
0050       double covGen = posDiff * posDiff;
0051       //       double covSym(covGen.LowerBlock());
0052       measCovar2 += weight * mixtureIter2->weight() * covGen;
0053     }
0054   }
0055
0056   double measCovar;
0057   if (weightSum < DBL_MIN) {
0058     std::cout << "MultiGaussianStateCombiner1D:: New state has total weight of 0." << std::endl;
0059     //     meanMean = SingleState::Vector(size,0);
0060     meanMean = 0.;
0061     measCovar = 0.;
0062     weightSum = 0.;
0063   } else {
0064     weightSum = 1. / weightSum;
0065     meanMean *= weightSum;
0066     measCovar1 *= weightSum;
0067     measCovar2 *= weightSum * weightSum;
0068     measCovar = measCovar1 + measCovar2;
0069     //     measCovar = measCovar1/weightSum + measCovar2/weightSum/weightSum;
0070   }
0071
0072   return SingleGaussianState1D(meanMean, measCovar, weightSum);
0073 }