KinematicFit/src/VertexKinematicConstraintT.cc

0001 #include "RecoVertex/KinematicFit/interface/VertexKinematicConstraintT.h"
0002 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 VertexKinematicConstraintT::VertexKinematicConstraintT() {}
0006
0007 VertexKinematicConstraintT::~VertexKinematicConstraintT() {}
0008
0009 void VertexKinematicConstraintT::init(const std::vector<KinematicState>& states,
0010                                       const GlobalPoint& ipoint,
0011                                       const GlobalVector& fieldValue) {
0012   int num = states.size();
0013   if (num != 2)
0014     throw VertexException("VertexKinematicConstraintT !=2 states passed");
0015   double mfz = fieldValue.z();
0016
0017   int j = 0;
0018   for (std::vector<KinematicState>::const_iterator i = states.begin(); i != states.end(); i++) {
0019     mom[j] = i->globalMomentum();
0020     dpos[j] = ipoint - i->globalPosition();
0021     a_i[j] = -i->particleCharge() * mfz;
0022
0023     double pvx = mom[j].x() - a_i[j] * dpos[j].y();
0024     double pvy = mom[j].y() + a_i[j] * dpos[j].x();
0025     double pvt2 = pvx * pvx + pvy * pvy;
0026     novera[j] = (dpos[j].x() * mom[j].x() + dpos[j].y() * mom[j].y());
0027     n[j] = a_i[j] * novera[j];
0028     m[j] = (pvx * mom[j].x() + pvy * mom[j].y());
0029     k[j] = -mom[j].z() / (mom[j].perp2() * pvt2);
0030     delta[j] = std::atan2(n[j], m[j]);
0031
0032     ++j;
0033   }
0034 }
0035
0036 void VertexKinematicConstraintT::fillValue() const {
0037   //it is 2 equations per track
0038   for (int j = 0; j != 2; ++j) {
0039     if (a_i[j] != 0) {
0040       //vector of values
0041       super::vl(j * 2) = dpos[j].y() * mom[j].x() - dpos[j].x() * mom[j].y() -
0042                          a_i[j] * (dpos[j].x() * dpos[j].x() + dpos[j].y() * dpos[j].y()) * 0.5;
0043       super::vl(j * 2 + 1) = dpos[j].z() - mom[j].z() * delta[j] / a_i[j];
0044     } else {
0045       //neutral particle
0046       double pt2Inverse = 1. / mom[j].perp2();
0047       super::vl(j * 2) = dpos[j].y() * mom[j].x() - dpos[j].x() * mom[j].y();
0048       super::vl(j * 2 + 1) =
0049           dpos[j].z() - mom[j].z() * ((dpos[j].x() * mom[j].x() + dpos[j].y() * mom[j].y()) * pt2Inverse);
0050     }
0051   }
0052 }
0053
0054 void VertexKinematicConstraintT::fillParametersDerivative() const {
0055   ROOT::Math::SMatrix<double, 2, 7> el_part_d;
0056   for (int j = 0; j != 2; ++j) {
0057     if (a_i[j] != 0) {
0058       //charged particle
0059
0060       //D Jacobian matrix
0061       el_part_d(0, 0) = mom[j].y() + a_i[j] * dpos[j].x();
0062       el_part_d(0, 1) = -mom[j].x() + a_i[j] * dpos[j].y();
0063       el_part_d(1, 0) = -k[j] * (m[j] * mom[j].x() - n[j] * mom[j].y());
0064       el_part_d(1, 1) = -k[j] * (m[j] * mom[j].y() + n[j] * mom[j].x());
0065       el_part_d(1, 2) = -1.;
0066       el_part_d(0, 3) = dpos[j].y();
0067       el_part_d(0, 4) = -dpos[j].x();
0068       el_part_d(1, 3) = k[j] * (m[j] * dpos[j].x() - novera[j] * (2 * mom[j].x() - a_i[j] * dpos[j].y()));
0069       el_part_d(1, 4) = k[j] * (m[j] * dpos[j].y() - novera[j] * (2 * mom[j].y() + a_i[j] * dpos[j].x()));
0070       el_part_d(1, 5) = -delta[j] / a_i[j];
0071       super::jac_d().Place_at(el_part_d, j * 2, j * 7);
0072     } else {
0073       //neutral particle
0074       double pt2Inverse = 1. / mom[j].perp2();
0075       el_part_d(0, 0) = mom[j].y();
0076       el_part_d(0, 1) = -mom[j].x();
0077       el_part_d(1, 0) = mom[j].x() * (mom[j].z() * pt2Inverse);
0078       el_part_d(1, 1) = mom[j].y() * (mom[j].z() * pt2Inverse);
0079       el_part_d(1, 2) = -1.;
0080       el_part_d(0, 3) = dpos[j].y();
0081       el_part_d(0, 4) = -dpos[j].x();
0082       el_part_d(1, 3) = 2 * (dpos[j].x() * mom[j].x() + dpos[j].y() * mom[j].y()) * pt2Inverse * mom[j].x() *
0083                             (mom[j].z() * pt2Inverse) -
0084                         dpos[j].x() * (mom[j].z() * pt2Inverse);
0085       el_part_d(1, 4) = 2 * (dpos[j].x() * mom[j].x() + dpos[j].y() * mom[j].y()) * pt2Inverse * mom[j].y() *
0086                             (mom[j].z() * pt2Inverse) -
0087                         dpos[j].x() * (mom[j].z() * pt2Inverse);
0088       el_part_d(1, 5) = -(dpos[j].x() * mom[j].x() + dpos[j].y() * mom[j].y()) * pt2Inverse;
0089       super::jac_d().Place_at(el_part_d, j * 2, j * 7);
0090     }
0091   }
0092 }
0093
0094 void VertexKinematicConstraintT::fillPositionDerivative() const {
0095   ROOT::Math::SMatrix<double, 2, 3> el_part_e;
0096   for (int j = 0; j != 2; ++j) {
0097     if (a_i[j] != 0) {
0098       //charged particle
0099
0100       //E jacobian matrix
0101       el_part_e(0, 0) = -(mom[j].y() + a_i[j] * dpos[j].x());
0102       el_part_e(0, 1) = mom[j].x() - a_i[j] * dpos[j].y();
0103       el_part_e(1, 0) = k[j] * (m[j] * mom[j].x() - n[j] * mom[j].y());
0104       el_part_e(1, 1) = k[j] * (m[j] * mom[j].y() + n[j] * mom[j].x());
0105       el_part_e(1, 2) = 1;
0106       super::jac_e().Place_at(el_part_e, 2 * j, 0);
0107     } else {
0108       //neutral particle
0109       double pt2Inverse = 1. / mom[j].perp2();
0110       el_part_e(0, 0) = -mom[j].y();
0111       el_part_e(0, 1) = mom[j].x();
0112       el_part_e(1, 0) = -mom[j].x() * mom[j].z() * pt2Inverse;
0113       el_part_e(1, 1) = -mom[j].y() * mom[j].z() * pt2Inverse;
0114       el_part_e(1, 2) = 1;
0115       super::jac_e().Place_at(el_part_e, 2 * j, 0);
0116     }
0117   }
0118 }
0119
0120 int VertexKinematicConstraintT::numberOfEquations() const { return 2; }