GsfTools/interface/MultiGaussianStateAssembler.icc

0001 #include "TrackingTools/GsfTools/interface/GaussianStateLessWeight.h"
0002 #include "FWCore/Utilities/interface/Exception.h"
0003
0004 template <unsigned int N>
0005 void MultiGaussianStateAssembler<N>::addState(const SingleStatePtr& sgs) {
0006   //
0007   // refuse to add states after combination has been done
0008   //
0009   if (combinationDone)
0010     throw cms::Exception("LogicError") << "MultiGaussianStateAssembler: trying to add states after combination";
0011
0012   theValidWeightSum += sgs->weight();
0013
0014   theStates.push_back(sgs);
0015 }
0016
0017 template <unsigned int N>
0018 void MultiGaussianStateAssembler<N>::addState(const MultiState& mgs) {
0019   //
0020   // refuse to add states after combination has been done
0021   //
0022   if (combinationDone)
0023     throw cms::Exception("LogicError") << "MultiGaussianStateAssembler: trying to add states after combination";
0024   //
0025   // Add components (i.e. state to be added can be single or multi state)
0026   //
0027   SingleStateContainer components(mgs.components());
0028   addStateVector(components);
0029 }
0030
0031 template <unsigned int N>
0032 void MultiGaussianStateAssembler<N>::addStateVector(const SingleStateContainer& states) {
0033   //
0034   // refuse to add states after combination has been done
0035   //
0036   if (combinationDone)
0037     throw cms::Exception("LogicError") << "MultiGaussianStateAssembler: trying to add states after combination";
0038   //
0039   // sum up weights (all components are supposed to be valid!!!)
0040   //
0041   double sum = 0.;
0042   for (auto const& is : states)
0043     sum += (*is).weight();
0044   theValidWeightSum += sum;
0045   //
0046   // add to vector of states
0047   //
0048   theStates.insert(theStates.end(), states.begin(), states.end());
0049 }
0050
0051 template <unsigned int N>
0052 MultiGaussianState<N> MultiGaussianStateAssembler<N>::combinedState() {
0053   //
0054   // Prepare resulting state vector
0055   prepareCombinedState();
0056   return MultiState(theStates);
0057 }
0058
0059 template <unsigned int N>
0060 MultiGaussianState<N> MultiGaussianStateAssembler<N>::combinedState(const float newWeight) {
0061   //
0062   // Prepare resulting state vector
0063   //
0064   prepareCombinedState();
0065   //
0066   // return reweighted state
0067   //
0068   return reweightedCombinedState(newWeight);
0069 }
0070
0071 template <unsigned int N>
0072 bool MultiGaussianStateAssembler<N>::prepareCombinedState() {
0073   //
0074   // remaining part to be done only once
0075   //
0076   if (combinationDone)
0077     return true;
0078   combinationDone = true;
0079   //
0080   // Remove states with negligible weights
0081   //
0082   removeSmallWeights();
0083   return !theStates.empty();
0084 }
0085
0086 template <unsigned int N>
0087 MultiGaussianState<N> MultiGaussianStateAssembler<N>::reweightedCombinedState(const double newWeight) const {
0088   //
0089   // scaling factor
0090   //
0091   double factor = theValidWeightSum > 0. ? newWeight / theValidWeightSum : 1;
0092   //
0093   // create new vector of states & combined state
0094   //
0095   SingleStateContainer reweightedStates;
0096   reweightedStates.reserve(theStates.size());
0097   for (auto const& is : theStates) {
0098     auto oldWeight = (*is).weight();
0099     reweightedStates.emplace_back((*is).mean(), (*is).covariance(), factor * oldWeight);
0100   }
0101   return MultiState(reweightedStates);
0102 }
0103
0104 template <unsigned int N>
0105 void MultiGaussianStateAssembler<N>::removeSmallWeights() {
0106   //
0107   // check total weight
0108   //
0109   if (theValidWeightSum == 0.) {
0110     theStates.clear();
0111     return;
0112   }
0113   theStates.erase(std::remove_if(theStates.begin(),
0114                                  theStates.end(),
0115                                  [&](typename SingleStateContainer::value_type const& s) {
0116                                    return s->weight() < minFractionalWeight * theValidWeightSum;
0117                                  }),
0118                   theStates.end());
0119 }