File indexing completed on 2024-04-06 12:29:12
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
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
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
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
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
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
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
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
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
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
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; }