KinematicFit/src/VertexKinematicConstraint.cc

0001 #include "RecoVertex/KinematicFit/interface/VertexKinematicConstraint.h"
0002 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 VertexKinematicConstraint::VertexKinematicConstraint() {}
0006
0007 VertexKinematicConstraint::~VertexKinematicConstraint() {}
0008
0009 AlgebraicVector VertexKinematicConstraint::value(const std::vector<KinematicState>& states,
0010                                                  const GlobalPoint& point) const {
0011   int num = states.size();
0012   if (num < 2)
0013     throw VertexException("VertexKinematicConstraint::<2 states passed");
0014
0015   //it is 2 equations per track
0016   AlgebraicVector vl(2 * num, 0);
0017   int num_r = 0;
0018   for (std::vector<KinematicState>::const_iterator i = states.begin(); i != states.end(); i++) {
0019     TrackCharge ch = i->particleCharge();
0020     GlobalVector mom = i->globalMomentum();
0021     GlobalPoint pos = i->globalPosition();
0022     double d_x = point.x() - pos.x();
0023     double d_y = point.y() - pos.y();
0024     double d_z = point.z() - pos.z();
0025     double pt = mom.transverse();
0026
0027     if (ch != 0) {
0028       //charged particle
0029       double a_i = -ch * i->magneticField()->inInverseGeV(pos).z();
0030
0031       double pvx = mom.x() - a_i * d_y;
0032       double pvy = mom.y() + a_i * d_x;
0033       double n = a_i * (d_x * mom.x() + d_y * mom.y());
0034       double m = (pvx * mom.x() + pvy * mom.y());
0035       double delta = atan2(n, m);
0036
0037       //vector of values
0038       vl(num_r * 2 + 1) = d_y * mom.x() - d_x * mom.y() - a_i * (d_x * d_x + d_y * d_y) / 2;
0039       vl(num_r * 2 + 2) = d_z - mom.z() * delta / a_i;
0040     } else {
0041       //neutral particle
0042       vl(num_r * 2 + 1) = d_y * mom.x() - d_x * mom.y();
0043       vl(num_r * 2 + 2) = d_z - mom.z() * (d_x * mom.x() + d_y * mom.y()) / (pt * pt);
0044     }
0045     num_r++;
0046   }
0047   return vl;
0048 }
0049
0050 AlgebraicMatrix VertexKinematicConstraint::parametersDerivative(const std::vector<KinematicState>& states,
0051                                                                 const GlobalPoint& point) const {
0052   int num = states.size();
0053   if (num < 2)
0054     throw VertexException("VertexKinematicConstraint::<2 states passed");
0055   AlgebraicMatrix jac_d(2 * num, 7 * num);
0056   int num_r = 0;
0057   for (std::vector<KinematicState>::const_iterator i = states.begin(); i != states.end(); i++) {
0058     AlgebraicMatrix el_part_d(2, 7, 0);
0059     TrackCharge ch = i->particleCharge();
0060     GlobalVector mom = i->globalMomentum();
0061     GlobalPoint pos = i->globalPosition();
0062     double d_x = point.x() - pos.x();
0063     double d_y = point.y() - pos.y();
0064     double pt = mom.transverse();
0065
0066     if (ch != 0) {
0067       //charged particle
0068       double a_i = -ch * i->magneticField()->inInverseGeV(pos).z();
0069
0070       double pvx = mom.x() - a_i * d_y;
0071       double pvy = mom.y() + a_i * d_x;
0072       double pvt = sqrt(pvx * pvx + pvy * pvy);
0073       double novera = (d_x * mom.x() + d_y * mom.y());
0074       double n = a_i * novera;
0075       double m = (pvx * mom.x() + pvy * mom.y());
0076       double k = -mom.z() / (pvt * pvt * pt * pt);
0077       double delta = atan2(n, m);
0078
0079       //D Jacobian matrix
0080       el_part_d(1, 1) = mom.y() + a_i * d_x;
0081       el_part_d(1, 2) = -mom.x() + a_i * d_y;
0082       el_part_d(2, 1) = -k * (m * mom.x() - n * mom.y());
0083       el_part_d(2, 2) = -k * (m * mom.y() + n * mom.x());
0084       el_part_d(2, 3) = -1.;
0085       el_part_d(1, 4) = d_y;
0086       el_part_d(1, 5) = -d_x;
0087       el_part_d(2, 4) = k * (m * d_x - novera * (2 * mom.x() - a_i * d_y));
0088       el_part_d(2, 5) = k * (m * d_y - novera * (2 * mom.y() + a_i * d_x));
0089       el_part_d(2, 6) = -delta / a_i;
0090       jac_d.sub(num_r * 2 + 1, num_r * 7 + 1, el_part_d);
0091     } else {
0092       //neutral particle
0093       el_part_d(1, 1) = mom.y();
0094       el_part_d(1, 2) = -mom.x();
0095       el_part_d(2, 1) = mom.x() * mom.z() / (pt * pt);
0096       el_part_d(2, 2) = mom.y() * mom.z() / (pt * pt);
0097       el_part_d(2, 3) = -1.;
0098       el_part_d(1, 4) = d_y;
0099       el_part_d(1, 5) = -d_x;
0100       el_part_d(2, 4) =
0101           2 * (d_x * mom.x() + d_y * mom.y()) * mom.x() * mom.z() / (pt * pt * pt * pt) - mom.z() * d_x / (pt * pt);
0102       el_part_d(2, 5) =
0103           2 * (d_x * mom.x() + d_y * mom.y()) * mom.y() * mom.z() / (pt * pt * pt * pt) - mom.z() * d_y / (pt * pt);
0104       el_part_d(2, 6) = -(d_x * mom.x() + d_y * mom.y()) / (pt * pt);
0105       jac_d.sub(num_r * 2 + 1, num_r * 7 + 1, el_part_d);
0106     }
0107     num_r++;
0108   }
0109   return jac_d;
0110 }
0111
0112 AlgebraicMatrix VertexKinematicConstraint::positionDerivative(const std::vector<KinematicState>& states,
0113                                                               const GlobalPoint& point) const {
0114   int num = states.size();
0115   if (num < 2)
0116     throw VertexException("VertexKinematicConstraint::<2 states passed");
0117   AlgebraicMatrix jac_e(2 * num, 3);
0118   int num_r = 0;
0119   for (std::vector<KinematicState>::const_iterator i = states.begin(); i != states.end(); i++) {
0120     AlgebraicMatrix el_part_e(2, 3, 0);
0121     TrackCharge ch = i->particleCharge();
0122     GlobalVector mom = i->globalMomentum();
0123     GlobalPoint pos = i->globalPosition();
0124     double d_x = point.x() - pos.x();
0125     double d_y = point.y() - pos.y();
0126     double pt = mom.transverse();
0127
0128     if (ch != 0) {
0129       //charged particle
0130       double a_i = -ch * i->magneticField()->inInverseGeV(pos).z();
0131
0132       double pvx = mom.x() - a_i * d_y;
0133       double pvy = mom.y() + a_i * d_x;
0134       double pvt = sqrt(pvx * pvx + pvy * pvy);
0135       double n = a_i * (d_x * mom.x() + d_y * mom.y());
0136       double m = (pvx * mom.x() + pvy * mom.y());
0137       double k = -mom.z() / (pvt * pvt * pt * pt);
0138
0139       //E jacobian matrix
0140       el_part_e(1, 1) = -(mom.y() + a_i * d_x);
0141       el_part_e(1, 2) = mom.x() - a_i * d_y;
0142       el_part_e(2, 1) = k * (m * mom.x() - n * mom.y());
0143       el_part_e(2, 2) = k * (m * mom.y() + n * mom.x());
0144       el_part_e(2, 3) = 1;
0145       jac_e.sub(2 * num_r + 1, 1, el_part_e);
0146     } else {
0147       //neutral particle
0148       el_part_e(1, 1) = -mom.y();
0149       el_part_e(1, 2) = mom.x();
0150       el_part_e(2, 1) = -mom.x() * mom.z() / (pt * pt);
0151       el_part_e(2, 2) = -mom.y() * mom.z() / (pt * pt);
0152       el_part_e(2, 3) = 1;
0153       jac_e.sub(2 * num_r + 1, 1, el_part_e);
0154     }
0155     num_r++;
0156   }
0157   return jac_e;
0158 }
0159
0160 int VertexKinematicConstraint::numberOfEquations() const { return 2; }