GaussianSumVertexFit/src/GsfVertexWeightCalculator.cc

0001 #include "RecoVertex/GaussianSumVertexFit/interface/GsfVertexWeightCalculator.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "FWCore/Utilities/interface/isFinite.h"
0004
0005 #include <cfloat>
0006
0007 double GsfVertexWeightCalculator::calculate(const VertexState& oldVertex,
0008                                             const RefCountedLinearizedTrackState track,
0009                                             double cov) const {
0010   double previousWeight = oldVertex.weightInMixture() * track->weightInMixture();
0011   // Jacobians
0012   const AlgebraicMatrixN3& a = track->positionJacobian();
0013   const AlgebraicMatrixNM& b = track->momentumJacobian();
0014
0015   //track information
0016   AlgebraicVectorN trackParameters = track->predictedStateParameters();
0017   AlgebraicSymMatrixNN trackParametersError = track->predictedStateError();
0018   //vertex information
0019   AlgebraicSymMatrix33 oldVertexError = oldVertex.error().matrix();
0020   //Vertex position
0021   GlobalPoint oldVertexPosition = oldVertex.position();
0022   AlgebraicVector3 oldVertexCoord;
0023   oldVertexCoord[0] = oldVertexPosition.x();
0024   oldVertexCoord[1] = oldVertexPosition.y();
0025   oldVertexCoord[2] = oldVertexPosition.z();
0026
0027   // prior momentum information
0028   AlgebraicVector3 priorMomentum = track->predictedStateMomentumParameters();
0029
0030   AlgebraicSymMatrix33 priorMomentumCov;
0031   priorMomentumCov(0, 0) = 1.0E-9;
0032   priorMomentumCov(1, 1) = 1.0E-6;
0033   priorMomentumCov(2, 2) = 1.0E-6;
0034   priorMomentumCov *= cov;
0035
0036   AlgebraicVectorN diff = trackParameters - track->constantTerm() - a * oldVertexCoord - b * priorMomentum;
0037   track->checkParameters(diff);
0038   AlgebraicSymMatrixNN sigmaM =
0039       trackParametersError + ROOT::Math::Similarity(a, oldVertexError) + ROOT::Math::Similarity(b, priorMomentumCov);
0040
0041   double sigmaDet;
0042   sigmaM.Det(sigmaDet);
0043   int ifail = !sigmaM.Invert();
0044   if (ifail != 0) {
0045     edm::LogWarning("GsfVertexWeightCalculator") << "S matrix inversion failed";
0046     return -1.;
0047   }
0048
0049   double chi = ROOT::Math::Similarity(diff, sigmaM);
0050   ;  //SigmaM is now inverted !!!
0051   double weight = pow(2. * M_PI, -0.5 * 5) * sqrt(1. / sigmaDet) * exp(-0.5 * chi);
0052
0053   if (edm::isNotFinite(weight) || sigmaDet <= 0.) {
0054     edm::LogWarning("GsfVertexWeightCalculator") << "Weight is NaN";
0055     return -1.;
0056   }
0057
0058   return weight * previousWeight;
0059 }