GsfTools/interface/MultiGaussianStateCombiner.icc

0001 #include "FWCore/Utilities/interface/Exception.h"
0002
0003 #include <cfloat>
0004
0005 template <unsigned int N>
0006 typename MultiGaussianStateCombiner<N>::SingleStatePtr MultiGaussianStateCombiner<N>::combine(
0007     const MultiState& theState) const {
0008   return combine(theState.components());
0009 }
0010
0011 template <unsigned int N>
0012 typename MultiGaussianStateCombiner<N>::SingleStatePtr MultiGaussianStateCombiner<N>::combine(
0013     const VSC& theComponents) const {
0014   if (theComponents.empty()) {
0015     throw cms::Exception("LogicError") << "MultiGaussianStateCombiner:: state container to collapse is empty";
0016     return SingleStatePtr();
0017   }
0018
0019   const SingleStatePtr& firstState(theComponents.front());
0020   if (theComponents.size() == 1) {
0021     return firstState;
0022   }
0023
0024   //   int size = firstState.mean().num_row();
0025   typename SingleState::Vector meanMean;
0026   double weightSum = 0.;
0027   double weight;
0028   typename SingleState::Matrix measCovar1, measCovar2;
0029   for (auto mixtureIter1 = theComponents.begin(); mixtureIter1 != theComponents.end(); mixtureIter1++) {
0030     weight = (*mixtureIter1)->weight();
0031     weightSum += weight;
0032
0033     typename SingleState::Vector mean1 = (**mixtureIter1).mean();
0034     meanMean += weight * mean1;
0035     measCovar1 += weight * (**mixtureIter1).covariance();
0036
0037     for (auto mixtureIter2 = mixtureIter1 + 1; mixtureIter2 != theComponents.end(); mixtureIter2++) {
0038       typename SingleState::Vector posDiff = mean1 - (**mixtureIter2).mean();
0039       //
0040       // TensorProd yields a general matrix - need to convert to a symm. matrix
0041       ROOT::Math::SMatrix<double, N, N> covGen = ROOT::Math::TensorProd(posDiff, posDiff);
0042       typename SingleState::Matrix covSym(covGen.LowerBlock());
0043       measCovar2 += weight * (**mixtureIter2).weight() * covSym;
0044     }
0045   }
0046
0047   typename SingleState::Matrix measCovar;
0048   if (weightSum < DBL_MIN) {
0049     std::cout << "MultiGaussianStateCombiner:: New state has total weight of 0." << std::endl;
0050     //     meanMean = SingleState::Vector(size,0);
0051     meanMean *= 0.;
0052     weightSum = 0.;
0053   } else {
0054     auto wsInv = 1. / weightSum;
0055     meanMean *= wsInv;
0056     measCovar1 *= wsInv;
0057     measCovar2 *= wsInv * wsInv;
0058     measCovar = measCovar1 + measCovar2;
0059   }
0060
0061   return std::make_shared<SingleState>(meanMean, measCovar, weightSum);
0062 }