KinematicFit/src/MomentumKinematicConstraint.cc

0001 #include "RecoVertex/KinematicFit/interface/MomentumKinematicConstraint.h"
0002 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
0003
0004 MomentumKinematicConstraint::MomentumKinematicConstraint(const AlgebraicVector& momentum, const AlgebraicVector& dev) {
0005   if ((momentum.num_row() != 3) || (dev.num_row() != 3))
0006     throw VertexException("MomentumKinemticConstraint::Momentum or Deviation vector passed is not 3-dimensional");
0007   mm = momentum;
0008   AlgebraicVector dev_l(7, 0);
0009   dev_l(4) = dev(1) * dev(1);
0010   dev_l(5) = dev(2) * dev(2);
0011   dev_l(6) = dev(3) * dev(3);
0012   dd = dev_l;
0013 }
0014
0015 std::pair<AlgebraicVector, AlgebraicVector> MomentumKinematicConstraint::value(const AlgebraicVector& exPoint) const {
0016   if (exPoint.num_row() == 0)
0017     throw VertexException("MomentumKinematicConstraint::value requested for zero Linearization point");
0018
0019   //security check for extended cartesian parametrization
0020   int inSize = exPoint.num_row();
0021   if ((inSize % 7) != 0)
0022     throw VertexException("MomentumKinematicConstraint::linearization point has a wrong dimension");
0023   int nStates = inSize / 7;
0024   if (nStates != 1)
0025     throw VertexException("MomentumKinematicConstraint::Multistate refit is not foreseen for this constraint");
0026   AlgebraicVector pr = exPoint;
0027   AlgebraicVector vl(3, 0);
0028   vl(1) = pr(4) - mm(1);
0029   vl(2) = pr(5) - mm(2);
0030   vl(3) = pr(6) - mm(3);
0031   return std::pair<AlgebraicVector, AlgebraicVector>(vl, pr);
0032 }
0033
0034 std::pair<AlgebraicMatrix, AlgebraicVector> MomentumKinematicConstraint::derivative(
0035     const AlgebraicVector& exPoint) const {
0036   if (exPoint.num_row() == 0)
0037     throw VertexException("MomentumKinematicConstraint::derivative requested for zero Linearization point");
0038
0039   //security check for extended cartesian parametrization
0040   int inSize = exPoint.num_row();
0041   if ((inSize % 7) != 0)
0042     throw VertexException("MomentumKinematicConstraint::linearization point has a wrong dimension");
0043   int nStates = inSize / 7;
0044   if (nStates != 1)
0045     throw VertexException("MomentumKinematicConstraint::Multistate refit is not foreseen for this constraint");
0046
0047   AlgebraicVector pr = exPoint;
0048   AlgebraicMatrix dr(3, 7, 0);
0049   dr(1, 4) = 1.;
0050   dr(2, 5) = 1.;
0051   dr(3, 6) = 1.;
0052   return std::pair<AlgebraicMatrix, AlgebraicVector>(dr, pr);
0053 }
0054
0055 std::pair<AlgebraicVector, AlgebraicVector> MomentumKinematicConstraint::value(
0056     const std::vector<RefCountedKinematicParticle>& par) const {
0057   int nStates = par.size();
0058   if (nStates == 0)
0059     throw VertexException("MomentumKinematicConstraint::Empty vector of particles passed");
0060   if (nStates != 1)
0061     throw VertexException("MomentumKinematicConstraint::Multistate refit is not foreseen for this constraint");
0062   AlgebraicVector point = asHepVector<7>(par.front()->currentState().kinematicParameters().vector());
0063   AlgebraicVector vl(3, 0);
0064   vl(1) = point(4) - mm(1);
0065   vl(2) = point(5) - mm(2);
0066   vl(3) = point(6) - mm(3);
0067   return std::pair<AlgebraicVector, AlgebraicVector>(vl, point);
0068 }
0069
0070 std::pair<AlgebraicMatrix, AlgebraicVector> MomentumKinematicConstraint::derivative(
0071     const std::vector<RefCountedKinematicParticle>& par) const {
0072   int nStates = par.size();
0073   if (nStates == 0)
0074     throw VertexException("MomentumKinematicConstraint::Empty vector of particles passed");
0075   if (nStates != 1)
0076     throw VertexException("MomentumKinematicConstraint::Multistate refit is not foreseen for this constraint");
0077   AlgebraicVector point = asHepVector<7>(par.front()->currentState().kinematicParameters().vector());
0078   AlgebraicMatrix dr(3, 7, 0);
0079   dr(1, 4) = 1.;
0080   dr(2, 5) = 1.;
0081   dr(3, 6) = 1.;
0082   return std::pair<AlgebraicMatrix, AlgebraicVector>(dr, point);
0083 }
0084
0085 AlgebraicVector MomentumKinematicConstraint::deviations(int nStates) const {
0086   if (nStates == 0)
0087     throw VertexException("MomentumKinematicConstraint::Empty vector of particles passed");
0088   if (nStates != 1)
0089     throw VertexException("MomentumKinematicConstraint::Multistate refit is not foreseen for this constraint");
0090   AlgebraicVector res = dd;
0091   return res;
0092 }
0093
0094 int MomentumKinematicConstraint::numberOfEquations() const { return 3; }