CMS 3D CMS Logo

List of all members | Public Member Functions | Private Types | Private Attributes
PosteriorWeightsCalculator Class Reference

#include <PosteriorWeightsCalculator.h>

Public Member Functions

 PosteriorWeightsCalculator (const std::vector< TSOS > &mixture)
 
std::vector< double > weights (const TrackingRecHit &tsos) const
 Create random state. More...
 
template<unsigned int D>
std::vector< double > weights (const TrackingRecHit &tsos) const
 
 ~PosteriorWeightsCalculator ()
 

Private Types

typedef TrajectoryStateOnSurface TSOS
 

Private Attributes

std::vector< TSOSpredictedComponents
 

Detailed Description

Helper class which calculates the posterior weights of a Gaussian mixture given a prior (predicted) mixture and a RecHit. The prior is specified during construction time in the form of a vector of trajectory states.

Definition at line 13 of file PosteriorWeightsCalculator.h.

Member Typedef Documentation

◆ TSOS

Definition at line 15 of file PosteriorWeightsCalculator.h.

Constructor & Destructor Documentation

◆ PosteriorWeightsCalculator()

PosteriorWeightsCalculator::PosteriorWeightsCalculator ( const std::vector< TSOS > &  mixture)
inline

Definition at line 18 of file PosteriorWeightsCalculator.h.

18 : predictedComponents(mixture) {}

◆ ~PosteriorWeightsCalculator()

PosteriorWeightsCalculator::~PosteriorWeightsCalculator ( )
inline

Definition at line 20 of file PosteriorWeightsCalculator.h.

20 {}

Member Function Documentation

◆ weights() [1/2]

std::vector< double > PosteriorWeightsCalculator::weights ( const TrackingRecHit tsos) const

Create random state.

Definition at line 11 of file PosteriorWeightsCalculator.cc.

11  {
12  switch (recHit.dimension()) {
13  case 1:
14  return weights<1>(recHit);
15  case 2:
16  return weights<2>(recHit);
17  case 3:
18  return weights<3>(recHit);
19  case 4:
20  return weights<4>(recHit);
21  case 5:
22  return weights<5>(recHit);
23  }
24  throw cms::Exception("Error: rechit of size not 1,2,3,4,5");
25 }

References Exception, and rpcPointValidation_cfi::recHit.

Referenced by GsfChi2MeasurementEstimator::estimate(), GsfMultiStateUpdator::update(), and weights().

◆ weights() [2/2]

template<unsigned int D>
std::vector< double > PosteriorWeightsCalculator::weights ( const TrackingRecHit tsos) const

Definition at line 28 of file PosteriorWeightsCalculator.cc.

28  {
29  typedef typename AlgebraicROOTObject<D, 5>::Matrix MatD5;
30  typedef typename AlgebraicROOTObject<5, D>::Matrix Mat5D;
31  typedef typename AlgebraicROOTObject<D, D>::SymMatrix SMatDD;
32  typedef typename AlgebraicROOTObject<D>::Vector VecD;
33  using ROOT::Math::SMatrixNoInit;
34 
35  std::vector<double> weights;
36  if (predictedComponents.empty()) {
37  edm::LogError("EmptyPredictedComponents") << "a multi state is empty. cannot compute any weight.";
38  return weights;
39  }
40  weights.reserve(predictedComponents.size());
41 
42  std::vector<double> detRs;
43  detRs.reserve(predictedComponents.size());
44  std::vector<double> chi2s;
45  chi2s.reserve(predictedComponents.size());
46 
47  VecD r, rMeas;
48  SMatDD V(SMatrixNoInit{}), R(SMatrixNoInit{});
50  //
51  // calculate chi2 and determinant / component and find
52  // minimum / maximum of chi2
53  //
54  double chi2Min(DBL_MAX);
55  for (unsigned int i = 0; i < predictedComponents.size(); i++) {
56  KfComponentsHolder holder;
57  auto const& x = predictedComponents[i].localParameters().vector();
58  holder.template setup<D>(&r, &V, &p, &rMeas, &R, x, predictedComponents[i].localError().matrix());
59  recHit.getKfComponents(holder);
60 
61  r -= rMeas;
62  R += V;
63 
64  double detR;
65  if (!R.Det2(detR)) {
66  edm::LogError("PosteriorWeightsCalculator") << "PosteriorWeightsCalculator: determinant failed";
67  return std::vector<double>();
68  }
69  detRs.push_back(detR);
70 
71  bool ok = invertPosDefMatrix(R);
72  if (!ok) {
73  edm::LogError("PosteriorWeightsCalculator") << "PosteriorWeightsCalculator: inversion failed";
74  return std::vector<double>();
75  }
76  double chi2 = ROOT::Math::Similarity(r, R);
77  chi2s.push_back(chi2);
78  if (chi2 < chi2Min)
79  chi2Min = chi2;
80  }
81 
82  if (detRs.size() != predictedComponents.size() || chi2s.size() != predictedComponents.size()) {
83  edm::LogError("PosteriorWeightsCalculator") << "Problem in vector sizes";
84  return std::vector<double>();
85  }
86 
87  //
88  // calculate weights (extracting a common factor
89  // exp(-0.5*chi2Min) to avoid numerical problems
90  // during exponentation
91  //
92  double sumWeights(0.);
93  for (unsigned int i = 0; i < predictedComponents.size(); i++) {
94  double priorWeight = predictedComponents[i].weight();
95 
96  double chi2 = chi2s[i] - chi2Min;
97 
98  double tempWeight(0.);
99  if (detRs[i] > FLT_MIN) {
100  //
101  // Calculation of (non-normalised) weight. Common factors exp(-chi2Norm/2.) and
102  // 1./sqrt(2*pi*recHit.dimension()) have been omitted
103  //
104  tempWeight = priorWeight * std::sqrt(1. / detRs[i]) * std::exp(-0.5 * chi2);
105  } else {
106  LogDebug("GsfTrackFitters") << "PosteriorWeightsCalculator: detR < FLT_MIN !!";
107  }
108  weights.push_back(tempWeight);
109  sumWeights += tempWeight;
110  }
111 
112  if (sumWeights < DBL_MIN) {
113  LogDebug("GsfTrackFitters") << "PosteriorWeightsCalculator: sumWeight < DBL_MIN";
114  edm::LogError("PosteriorWeightsCalculator") << "PosteriorWeightsCalculator: sumWeight < DBL_MIN";
115  return std::vector<double>();
116  }
117 
118  if (weights.size() != predictedComponents.size()) {
119  edm::LogError("PosteriorWeightsCalculator") << "Problem in vector sizes (2)";
120  return std::vector<double>();
121  }
122  sumWeights = 1. / sumWeights;
123  for (auto& w : weights)
124  w *= sumWeights;
125  return weights;
126 }

References hltPixelTracks_cff::chi2, muonRecoAnalyzer_cfi::chi2Min, JetChargeProducer_cfi::exp, mps_fire::i, invertPosDefMatrix(), LogDebug, makeMuonMisalignmentScenario::matrix, convertSQLiteXML::ok, AlCaHLTBitMon_ParallelJobs::p, predictedComponents, dttmaxenums::R, alignCSCRings::r, rpcPointValidation_cfi::recHit, mathSSE::sqrt(), cms::cuda::V, w, weights(), and x.

Member Data Documentation

◆ predictedComponents

std::vector<TSOS> PosteriorWeightsCalculator::predictedComponents
private

Definition at line 27 of file PosteriorWeightsCalculator.h.

Referenced by weights().

mps_fire.i
i
Definition: mps_fire.py:428
makeMuonMisalignmentScenario.matrix
list matrix
Definition: makeMuonMisalignmentScenario.py:141
AlgebraicROOTObject::SymMatrix
ROOT::Math::SMatrix< double, D1, D1, ROOT::Math::MatRepSym< double, D1 > > SymMatrix
Definition: AlgebraicROOTObjects.h:68
muonRecoAnalyzer_cfi.chi2Min
chi2Min
Definition: muonRecoAnalyzer_cfi.py:34
DDAxes::x
hltPixelTracks_cff.chi2
chi2
Definition: hltPixelTracks_cff.py:25
convertSQLiteXML.ok
bool ok
Definition: convertSQLiteXML.py:98
rpcPointValidation_cfi.recHit
recHit
Definition: rpcPointValidation_cfi.py:7
w
const double w
Definition: UKUtility.cc:23
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
PosteriorWeightsCalculator::weights
std::vector< double > weights(const TrackingRecHit &tsos) const
Create random state.
Definition: PosteriorWeightsCalculator.cc:11
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:233
AlgebraicROOTObject::Matrix
ROOT::Math::SMatrix< double, D1, D2, ROOT::Math::MatRepStd< double, D1, D2 > > Matrix
Definition: AlgebraicROOTObjects.h:69
AlCaHLTBitMon_ParallelJobs.p
def p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
cms::cuda::V
uint32_t const T *__restrict__ const uint32_t *__restrict__ int32_t int Histo::index_type cudaStream_t V
Definition: HistoContainer.h:51
invertPosDefMatrix
bool invertPosDefMatrix(ROOT::Math::SMatrix< T, N, N, ROOT::Math::MatRepSym< T, N > > &m)
Definition: invertPosDefMatrix.h:10
edm::LogError
Log< level::Error, false > LogError
Definition: MessageLogger.h:123
KfComponentsHolder
Definition: KfComponentsHolder.h:13
PosteriorWeightsCalculator::predictedComponents
std::vector< TSOS > predictedComponents
Definition: PosteriorWeightsCalculator.h:27
alignCSCRings.r
r
Definition: alignCSCRings.py:93
Exception
Definition: hltDiff.cc:245
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
dttmaxenums::R
Definition: DTTMax.h:29
ProjectMatrix
Definition: ProjectMatrix.h:8
AlgebraicROOTObject::Vector
ROOT::Math::SVector< double, D1 > Vector
Definition: AlgebraicROOTObjects.h:67