TrackingTools/GsfTools/interface/MultiGaussianStateCombiner.icc - master - githubmirror/cmssw

TrackingTools/GsfTools/interface/MultiGaussianStateCombiner.icc

62 строки2 KB

Shahzad Malik Muzaffar

[RECONSTRUCTION] [LLVM14] Apply code checks

30 ноя 2022, 17:50

30 ноя 2022, 17:50aae3b64

#include "FWCore/Utilities/interface/Exception.h"

#include <cfloat>

template <unsigned int N>
typename MultiGaussianStateCombiner<N>::SingleStatePtr MultiGaussianStateCombiner<N>::combine(
    const MultiState& theState) const {
  return combine(theState.components());
}

template <unsigned int N>
typename MultiGaussianStateCombiner<N>::SingleStatePtr MultiGaussianStateCombiner<N>::combine(
    const VSC& theComponents) const {
  if (theComponents.empty()) {
    throw cms::Exception("LogicError") << "MultiGaussianStateCombiner:: state container to collapse is empty";
    return SingleStatePtr();
  }

  const SingleStatePtr& firstState(theComponents.front());
  if (theComponents.size() == 1) {
    return firstState;
  }

  //   int size = firstState.mean().num_row();
  typename SingleState::Vector meanMean;
  double weightSum = 0.;
  double weight;
  typename SingleState::Matrix measCovar1, measCovar2;
  for (auto mixtureIter1 = theComponents.begin(); mixtureIter1 != theComponents.end(); mixtureIter1++) {
    weight = (*mixtureIter1)->weight();
    weightSum += weight;

    typename SingleState::Vector mean1 = (**mixtureIter1).mean();
    meanMean += weight * mean1;
    measCovar1 += weight * (**mixtureIter1).covariance();

    for (auto mixtureIter2 = mixtureIter1 + 1; mixtureIter2 != theComponents.end(); mixtureIter2++) {
      typename SingleState::Vector posDiff = mean1 - (**mixtureIter2).mean();
      //
      // TensorProd yields a general matrix - need to convert to a symm. matrix
      ROOT::Math::SMatrix<double, N, N> covGen = ROOT::Math::TensorProd(posDiff, posDiff);
      typename SingleState::Matrix covSym(covGen.LowerBlock());
      measCovar2 += weight * (**mixtureIter2).weight() * covSym;
    }
  }

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

  return std::make_shared<SingleState>(meanMean, measCovar, weightSum);
}