KinematicFit/src/FourMomentumKinematicConstraint.cc

0001 #include "RecoVertex/KinematicFit/interface/FourMomentumKinematicConstraint.h"
0002 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
0003
0004 FourMomentumKinematicConstraint::FourMomentumKinematicConstraint(const AlgebraicVector& momentum,
0005                                                                  const AlgebraicVector& deviation) {
0006   if ((momentum.num_row() != 4) || (deviation.num_row() != 4))
0007     throw VertexException("FourMomentumKinematicConstraint::vector of wrong size passed as 4-Momentum or Deviations");
0008   mm = momentum;
0009   AlgebraicVector deviation_l(7, 0);
0010   deviation_l(6) = momentum(3);
0011   deviation_l(5) = momentum(2);
0012   deviation_l(4) = momentum(1);
0013
0014   double mass_sq =
0015       momentum(4) * momentum(4) - momentum(3) * momentum(3) - momentum(2) * momentum(2) - momentum(1) * momentum(1);
0016
0017   if (mass_sq == 0.)
0018     throw VertexException("FourMomentumKinematicConstraint::the constraint vector passed corresponds to 0 mass");
0019   //deviation for mass calculated from deviations
0020   //of momentum
0021   deviation_l(7) =
0022       momentum(4) * momentum(4) * deviation(4) / mass_sq + momentum(3) * momentum(3) * deviation(3) / mass_sq +
0023       momentum(2) * momentum(2) * deviation(2) / mass_sq + momentum(1) * momentum(1) * deviation(1) / mass_sq;
0024   //mass sigma because of the energy
0025
0026   dd = deviation_l;
0027 }
0028
0029 std::pair<AlgebraicVector, AlgebraicVector> FourMomentumKinematicConstraint::value(
0030     const AlgebraicVector& exPoint) const {
0031   if (exPoint.num_row() == 0)
0032     throw VertexException("FourMomentumKinematicConstraint::value requested for zero Linearization point");
0033
0034   //security check for extended cartesian parametrization
0035   int inSize = exPoint.num_row();
0036   if ((inSize % 7) != 0)
0037     throw VertexException("FourMomentumKinematicConstraint::linearization point has a wrong dimension");
0038   int nStates = inSize / 7;
0039   if (nStates != 1)
0040     throw VertexException("FourMomentumKinematicConstraint::Multi particle refit is not supported in this version");
0041   AlgebraicVector pr = exPoint;
0042   AlgebraicVector vl(4, 0);
0043   vl(1) += (pr(4) - mm(1));
0044   vl(2) += (pr(5) - mm(2));
0045   vl(3) += (pr(6) - mm(3));
0046   vl(7) += (sqrt(pr(4) * pr(4) + pr(5) * pr(5) + pr(6) * pr(6) + pr(7) * pr(7)) - mm(4));
0047
0048   return std::pair<AlgebraicVector, AlgebraicVector>(vl, pr);
0049 }
0050
0051 std::pair<AlgebraicMatrix, AlgebraicVector> FourMomentumKinematicConstraint::derivative(
0052     const AlgebraicVector& exPoint) const {
0053   if (exPoint.num_row() == 0)
0054     throw VertexException("FourMomentumKinematicConstraint::value requested for zero Linearization point");
0055
0056   //security check for extended cartesian parametrization
0057   int inSize = exPoint.num_row();
0058   if ((inSize % 7) != 0)
0059     throw VertexException("FourMomentumKinematicConstraint::linearization point has a wrong dimension");
0060   int nStates = inSize / 7;
0061   if (nStates != 1)
0062     throw VertexException("FourMomentumKinematicConstraint::Multi particle refit is not supported in this version");
0063   AlgebraicVector pr = exPoint;
0064   AlgebraicMatrix dr(4, 7, 0);
0065
0066   dr(1, 4) = 1.;
0067   dr(2, 5) = 1.;
0068   dr(3, 6) = 1.;
0069   dr(4, 7) = pr(7) / sqrt(pr(4) * pr(4) + pr(5) * pr(5) + pr(6) * pr(6) + pr(7) * pr(7));
0070
0071   return std::pair<AlgebraicMatrix, AlgebraicVector>(dr, pr);
0072 }
0073
0074 std::pair<AlgebraicVector, AlgebraicVector> FourMomentumKinematicConstraint::value(
0075     const std::vector<RefCountedKinematicParticle>& par) const {
0076   int nStates = par.size();
0077   if (nStates == 0)
0078     throw VertexException("FourMomentumKinematicConstraint::Empty vector of particles passed");
0079   if (nStates != 1)
0080     throw VertexException("FourMomentumKinematicConstraint::Multi particle refit is not supported in this version");
0081   AlgebraicVector pr = asHepVector<7>(par.front()->currentState().kinematicParameters().vector());
0082   AlgebraicVector vl(4, 0);
0083
0084   vl(1) = pr(4) - mm(1);
0085   vl(2) = pr(5) - mm(2);
0086   vl(3) = pr(6) - mm(3);
0087   vl(7) = sqrt(pr(4) * pr(4) + pr(5) * pr(5) + pr(6) * pr(6) + pr(7) * pr(7)) - mm(4);
0088
0089   return std::pair<AlgebraicVector, AlgebraicVector>(vl, pr);
0090 }
0091
0092 std::pair<AlgebraicMatrix, AlgebraicVector> FourMomentumKinematicConstraint::derivative(
0093     const std::vector<RefCountedKinematicParticle>& par) const {
0094   int nStates = par.size();
0095   if (nStates == 0)
0096     throw VertexException("FourMomentumKinematicConstraint::Empty vector of particles passed");
0097   if (nStates != 1)
0098     throw VertexException("FourMomentumKinematicConstraint::Multi particle refit is not supported in this version");
0099   AlgebraicMatrix dr(4, 7, 0);
0100
0101   AlgebraicVector pr = asHepVector<7>(par.front()->currentState().kinematicParameters().vector());
0102   dr(1, 4) = 1.;
0103   dr(2, 5) = 1.;
0104   dr(3, 6) = 1.;
0105   dr(4, 7) = -pr(7) / sqrt(pr(4) * pr(4) + pr(5) * pr(5) + pr(6) * pr(6) + pr(7) * pr(7));
0106
0107   return std::pair<AlgebraicMatrix, AlgebraicVector>(dr, pr);
0108 }
0109
0110 AlgebraicVector FourMomentumKinematicConstraint::deviations(int nStates) const {
0111   if (nStates == 0)
0112     throw VertexException("FourMomentumKinematicConstraint::Empty vector of particles passed");
0113   if (nStates != 1)
0114     throw VertexException("FourMomentumKinematicConstraint::Multi particle refit is not supported in this version");
0115   AlgebraicVector res = dd;
0116   return res;
0117 }
0118
0119 int FourMomentumKinematicConstraint::numberOfEquations() const { return 4; }