KinematicFit/src/TwoTrackMassKinematicConstraint.cc

0001 #include "RecoVertex/KinematicFit/interface/TwoTrackMassKinematicConstraint.h"
0002 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
0003
0004 AlgebraicVector TwoTrackMassKinematicConstraint::value(const std::vector<KinematicState>& states,
0005                                                        const GlobalPoint& point) const {
0006   if (states.size() < 2)
0007     throw VertexException("TwoTrackMassKinematicConstraint::<2 states passed");
0008
0009   AlgebraicVector res(1, 0);
0010
0011   double a_1 = -states[0].particleCharge() * states[0].magneticField()->inInverseGeV(states[0].globalPosition()).z();
0012   double a_2 = -states[1].particleCharge() * states[1].magneticField()->inInverseGeV(states[1].globalPosition()).z();
0013
0014   AlgebraicVector7 p1 = states[0].kinematicParameters().vector();
0015   AlgebraicVector7 p2 = states[1].kinematicParameters().vector();
0016
0017   double p1vx = p1(3) - a_1 * (point.y() - p1(1));
0018   double p1vy = p1(4) + a_1 * (point.x() - p1(0));
0019   double p1vz = p1(5);
0020   ParticleMass m1 = p1(6);
0021
0022   double p2vx = p2(3) - a_2 * (point.y() - p2(1));
0023   double p2vy = p2(4) + a_2 * (point.x() - p2(0));
0024   double p2vz = p2(5);
0025   ParticleMass m2 = p2(6);
0026
0027   double j_energy = sqrt(p1(3) * p1(3) + p1(4) * p1(4) + p1(5) * p1(5) + m1 * m1) +
0028                     sqrt(p2(3) * p2(3) + p2(4) * p2(4) + p2(5) * p2(5) + m2 * m2);
0029
0030   double j_m = (p1vx + p2vx) * (p1vx + p2vx) + (p1vy + p2vy) * (p1vy + p2vy) + (p1vz + p2vz) * (p1vz + p2vz);
0031
0032   res(1) = j_energy * j_energy - j_m - mass * mass;
0033   return res;
0034 }
0035
0036 AlgebraicMatrix TwoTrackMassKinematicConstraint::parametersDerivative(const std::vector<KinematicState>& states,
0037                                                                       const GlobalPoint& point) const {
0038   int n_st = states.size();
0039   if (n_st < 2)
0040     throw VertexException("TwoTrackMassKinematicConstraint::<2 states passed");
0041
0042   AlgebraicMatrix res(1, n_st * 7, 0);
0043
0044   double a_1 = -states[0].particleCharge() * states[0].magneticField()->inInverseGeV(states[0].globalPosition()).z();
0045   double a_2 = -states[1].particleCharge() * states[1].magneticField()->inInverseGeV(states[1].globalPosition()).z();
0046
0047   AlgebraicVector7 p1 = states[0].kinematicParameters().vector();
0048   AlgebraicVector7 p2 = states[1].kinematicParameters().vector();
0049
0050   double p1vx = p1(3) - a_1 * (point.y() - p1(1));
0051   double p1vy = p1(4) + a_1 * (point.x() - p1(0));
0052   double p1vz = p1(5);
0053   ParticleMass m1 = p1(6);
0054
0055   double p2vx = p2(3) - a_2 * (point.y() - p2(1));
0056   double p2vy = p2(4) + a_2 * (point.x() - p2(0));
0057   double p2vz = p2(5);
0058   ParticleMass m2 = p2(6);
0059
0060   double e1 = sqrt(p1(3) * p1(3) + p1(4) * p1(4) + p1(5) * p1(5) + m1 * m1);
0061   double e2 = sqrt(p2(3) * p2(3) + p2(4) * p2(4) + p2(5) * p2(5) + m2 * m2);
0062
0063   //x1 and x2 derivatives: 1st and 8th elements
0064   res(1, 1) = 2 * a_1 * (p2vy + p1vy);
0065   res(1, 8) = 2 * a_2 * (p2vy + p1vy);
0066
0067   //y1 and y2 derivatives: 2nd and 9th elements:
0068   res(1, 2) = -2 * a_1 * (p1vx + p2vx);
0069   res(1, 9) = -2 * a_2 * (p2vx + p1vx);
0070
0071   //z1 and z2 components: 3d and 10th elmnets stay 0:
0072   res(1, 3) = 0.;
0073   res(1, 10) = 0.;
0074
0075   //px1 and px2 components: 4th and 11th elements:
0076   res(1, 4) = 2 * (1 + e2 / e1) * p1(3) - 2 * (p1vx + p2vx);
0077   res(1, 11) = 2 * (1 + e1 / e2) * p2(3) - 2 * (p1vx + p2vx);
0078
0079   //py1 and py2 components: 5th and 12 elements:
0080   res(1, 5) = 2 * (1 + e2 / e1) * p1(4) - 2 * (p1vy + p2vy);
0081   res(1, 12) = 2 * (1 + e1 / e2) * p2(4) - 2 * (p2vy + p1vy);
0082
0083   //pz1 and pz2 components: 6th and 13 elements:
0084   res(1, 6) = 2 * (1 + e2 / e1) * p1(5) - 2 * (p1vz + p2vz);
0085   res(1, 13) = 2 * (1 + e1 / e2) * p2(5) - 2 * (p2vz + p1vz);
0086
0087   //mass components: 7th and 14th elements:
0088   res(1, 7) = 2 * m1 * (1 + e2 / e1);
0089   res(1, 14) = 2 * m2 * (1 + e1 / e2);
0090
0091   return res;
0092 }
0093
0094 AlgebraicMatrix TwoTrackMassKinematicConstraint::positionDerivative(const std::vector<KinematicState>& states,
0095                                                                     const GlobalPoint& point) const {
0096   AlgebraicMatrix res(1, 3, 0);
0097   if (states.size() < 2)
0098     throw VertexException("TwoTrackMassKinematicConstraint::<2 states passed");
0099
0100   double a_1 = -states[0].particleCharge() * states[0].magneticField()->inInverseGeV(states[0].globalPosition()).z();
0101   double a_2 = -states[1].particleCharge() * states[1].magneticField()->inInverseGeV(states[1].globalPosition()).z();
0102
0103   AlgebraicVector7 p1 = states[0].kinematicParameters().vector();
0104   AlgebraicVector7 p2 = states[1].kinematicParameters().vector();
0105
0106   double p1vx = p1(3) - a_1 * (point.y() - p1(1));
0107   double p1vy = p1(4) + a_1 * (point.x() - p1(0));
0108
0109   double p2vx = p2(3) - a_2 * (point.y() - p2(1));
0110   double p2vy = p2(4) + a_2 * (point.x() - p2(0));
0111
0112   //xv component
0113   res(1, 1) = -2 * (p1vy + p2vy) * (a_1 + a_2);
0114
0115   //yv component
0116   res(1, 2) = 2 * (p1vx + p2vx) * (a_1 + a_2);
0117
0118   //zv component
0119   res(1, 3) = 0.;
0120
0121   return res;
0122 }
0123
0124 int TwoTrackMassKinematicConstraint::numberOfEquations() const { return 1; }