00001 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexTrackUpdator.h"
00002 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
00003
00004 #include "RecoVertex/VertexPrimitives/interface/VertexTrack.h"
00005 #include "RecoVertex/VertexPrimitives/interface/CachingVertex.h"
00006 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
00007
00008
00009 template <unsigned int N>
00010 typename CachingVertex<N>::RefCountedVertexTrack
00011 KalmanVertexTrackUpdator<N>::update
00012 (const CachingVertex<N> & vertex , RefCountedVertexTrack track) const
00013
00014 {
00015 trackMatrixPair thePair =
00016 trackRefit(vertex.vertexState(), track->linearizedTrack(), track->weight() );
00017
00018 VertexState rVert = updator.positionUpdate (vertex.vertexState(), track->linearizedTrack(),
00019 track->weight(), -1);
00020
00021 std::pair<bool, double> result = helper.trackParameterChi2(track->linearizedTrack(), thePair.first);
00022 float smoothedChi2 = helper.vertexChi2(rVert, vertex.vertexState()) + result.second;
00023
00024 return theVTFactory.vertexTrack(track->linearizedTrack(),
00025 vertex.vertexState(), thePair.first, smoothedChi2, thePair.second,
00026 track->weight());
00027 }
00028
00029 template <unsigned int N>
00030 typename KalmanVertexTrackUpdator<N>::trackMatrixPair
00031 KalmanVertexTrackUpdator<N>::trackRefit(const VertexState & vertex,
00032 KalmanVertexTrackUpdator<N>::RefCountedLinearizedTrackState linTrackState,
00033 float weight) const
00034
00035 {
00036 typedef ROOT::Math::SVector<double,N> AlgebraicVectorN;
00037 typedef ROOT::Math::SVector<double,N-2> AlgebraicVectorM;
00038 typedef ROOT::Math::SMatrix<double,N,3,ROOT::Math::MatRepStd<double,N,3> > AlgebraicMatrixN3;
00039 typedef ROOT::Math::SMatrix<double,N,N-2,ROOT::Math::MatRepStd<double,N,N-2> > AlgebraicMatrixNM;
00040 typedef ROOT::Math::SMatrix<double,N-2,3,ROOT::Math::MatRepStd<double,N-2,3> > AlgebraicMatrixM3;
00041 typedef ROOT::Math::SMatrix<double,N,N,ROOT::Math::MatRepSym<double,N> > AlgebraicSymMatrixNN;
00042
00043 typedef ROOT::Math::SMatrix<double,N+1,N+1,ROOT::Math::MatRepStd<double,N+1,N+1> > AlgebraicMatrixOO;
00044 typedef ROOT::Math::SMatrix<double,N-2,N-2,ROOT::Math::MatRepSym<double,N-2> > AlgebraicSymMatrixMM;
00045
00046
00047 GlobalPoint vertexPosition = vertex.position();
00048
00049 AlgebraicVector3 vertexCoord;
00050 vertexCoord(0) = vertexPosition.x();
00051 vertexCoord(1) = vertexPosition.y();
00052 vertexCoord(2) = vertexPosition.z();
00053 AlgebraicSymMatrix33 vertexErrorMatrix = vertex.error().matrix_new();
00054
00055
00056 const AlgebraicMatrixN3 & a = linTrackState->positionJacobian();
00057 const AlgebraicMatrixNM & b = linTrackState->momentumJacobian();
00058
00059
00060
00061
00062 int ifail;
00063 AlgebraicSymMatrixNN trackParametersWeight =
00064 linTrackState->predictedStateWeight(ifail);
00065
00066 AlgebraicSymMatrixMM s = ROOT::Math::SimilarityT(b,trackParametersWeight);
00067
00068 ifail = ! s.Invert();
00069 if(ifail !=0) throw VertexException
00070 ("KalmanVertexTrackUpdator::S matrix inversion failed");
00071
00072 AlgebraicVectorM newTrackMomentumP = s * (ROOT::Math::Transpose(b)) * trackParametersWeight *
00073 (linTrackState->predictedStateParameters() - linTrackState->constantTerm() - a*vertexCoord);
00074
00075 AlgebraicMatrix3M refittedPositionMomentumConvariance =
00076 -vertexErrorMatrix * (ROOT::Math::Transpose(a)) * trackParametersWeight * b * s;
00077
00078 AlgebraicSymMatrixMM refittedMomentumConvariance = s/weight +
00079 ROOT::Math::SimilarityT(refittedPositionMomentumConvariance, vertex.weight().matrix_new());
00080
00081
00082
00083 AlgebraicMatrixOO covMatrix;
00084 covMatrix.Place_at(refittedPositionMomentumConvariance, 0, 3);
00085 covMatrix.Place_at(ROOT::Math::Transpose(refittedPositionMomentumConvariance), 3, 0);
00086 covMatrix.Place_at(vertexErrorMatrix, 0, 0);
00087 covMatrix.Place_at(refittedMomentumConvariance, 3 ,3);
00088
00089 AlgebraicSymMatrixOO covSymMatrix(covMatrix.LowerBlock());
00090
00091 RefCountedRefittedTrackState refittedTrackState = linTrackState->
00092 createRefittedTrackState(vertexPosition, newTrackMomentumP, covSymMatrix);
00093
00094 return trackMatrixPair(refittedTrackState, covSymMatrix);
00095
00096 }
00097
00098 template class KalmanVertexTrackUpdator<5>;
00099 template class KalmanVertexTrackUpdator<6>;