VertexKinematicConstraint.cc

Go to the documentation of this file.

#include "RecoVertex/KinematicFit/interface/VertexKinematicConstraint.h"
#include "RecoVertex/VertexPrimitives/interface/VertexException.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

VertexKinematicConstraint::VertexKinematicConstraint() {}

VertexKinematicConstraint::~VertexKinematicConstraint() {}

AlgebraicVector VertexKinematicConstraint::value(const std::vector<KinematicState>& states,
                                                 const GlobalPoint& point) const {
  int num = states.size();
  if (num < 2)
    throw VertexException("VertexKinematicConstraint::<2 states passed");

  //it is 2 equations per track
  AlgebraicVector vl(2 * num, 0);
  int num_r = 0;
  for (std::vector<KinematicState>::const_iterator i = states.begin(); i != states.end(); i++) {
    TrackCharge ch = i->particleCharge();
    GlobalVector mom = i->globalMomentum();
    GlobalPoint pos = i->globalPosition();
    double d_x = point.x() - pos.x();
    double d_y = point.y() - pos.y();
    double d_z = point.z() - pos.z();
    double pt = mom.transverse();

    if (ch != 0) {
      //charged particle
      double a_i = -ch * i->magneticField()->inInverseGeV(pos).z();

      double pvx = mom.x() - a_i * d_y;
      double pvy = mom.y() + a_i * d_x;
      double n = a_i * (d_x * mom.x() + d_y * mom.y());
      double m = (pvx * mom.x() + pvy * mom.y());
      double delta = atan2(n, m);

      //vector of values
      vl(num_r * 2 + 1) = d_y * mom.x() - d_x * mom.y() - a_i * (d_x * d_x + d_y * d_y) / 2;
      vl(num_r * 2 + 2) = d_z - mom.z() * delta / a_i;
    } else {
      //neutral particle
      vl(num_r * 2 + 1) = d_y * mom.x() - d_x * mom.y();
      vl(num_r * 2 + 2) = d_z - mom.z() * (d_x * mom.x() + d_y * mom.y()) / (pt * pt);
    }
    num_r++;
  }
  return vl;
}

AlgebraicMatrix VertexKinematicConstraint::parametersDerivative(const std::vector<KinematicState>& states,
                                                                const GlobalPoint& point) const {
  int num = states.size();
  if (num < 2)
    throw VertexException("VertexKinematicConstraint::<2 states passed");
  AlgebraicMatrix jac_d(2 * num, 7 * num);
  int num_r = 0;
  for (std::vector<KinematicState>::const_iterator i = states.begin(); i != states.end(); i++) {
    AlgebraicMatrix el_part_d(2, 7, 0);
    TrackCharge ch = i->particleCharge();
    GlobalVector mom = i->globalMomentum();
    GlobalPoint pos = i->globalPosition();
    double d_x = point.x() - pos.x();
    double d_y = point.y() - pos.y();
    double pt = mom.transverse();

    if (ch != 0) {
      //charged particle
      double a_i = -ch * i->magneticField()->inInverseGeV(pos).z();

      double pvx = mom.x() - a_i * d_y;
      double pvy = mom.y() + a_i * d_x;
      double pvt = sqrt(pvx * pvx + pvy * pvy);
      double novera = (d_x * mom.x() + d_y * mom.y());
      double n = a_i * novera;
      double m = (pvx * mom.x() + pvy * mom.y());
      double k = -mom.z() / (pvt * pvt * pt * pt);
      double delta = atan2(n, m);

      //D Jacobian matrix
      el_part_d(1, 1) = mom.y() + a_i * d_x;
      el_part_d(1, 2) = -mom.x() + a_i * d_y;
      el_part_d(2, 1) = -k * (m * mom.x() - n * mom.y());
      el_part_d(2, 2) = -k * (m * mom.y() + n * mom.x());
      el_part_d(2, 3) = -1.;
      el_part_d(1, 4) = d_y;
      el_part_d(1, 5) = -d_x;
      el_part_d(2, 4) = k * (m * d_x - novera * (2 * mom.x() - a_i * d_y));
      el_part_d(2, 5) = k * (m * d_y - novera * (2 * mom.y() + a_i * d_x));
      el_part_d(2, 6) = -delta / a_i;
      jac_d.sub(num_r * 2 + 1, num_r * 7 + 1, el_part_d);
    } else {
      //neutral particle
      el_part_d(1, 1) = mom.y();
      el_part_d(1, 2) = -mom.x();
      el_part_d(2, 1) = mom.x() * mom.z() / (pt * pt);
      el_part_d(2, 2) = mom.y() * mom.z() / (pt * pt);
      el_part_d(2, 3) = -1.;
      el_part_d(1, 4) = d_y;
      el_part_d(1, 5) = -d_x;
      el_part_d(2, 4) =
          2 * (d_x * mom.x() + d_y * mom.y()) * mom.x() * mom.z() / (pt * pt * pt * pt) - mom.z() * d_x / (pt * pt);
      el_part_d(2, 5) =
          2 * (d_x * mom.x() + d_y * mom.y()) * mom.y() * mom.z() / (pt * pt * pt * pt) - mom.z() * d_y / (pt * pt);
      el_part_d(2, 6) = -(d_x * mom.x() + d_y * mom.y()) / (pt * pt);
      jac_d.sub(num_r * 2 + 1, num_r * 7 + 1, el_part_d);
    }
    num_r++;
  }
  return jac_d;
}

AlgebraicMatrix VertexKinematicConstraint::positionDerivative(const std::vector<KinematicState>& states,
                                                              const GlobalPoint& point) const {
  int num = states.size();
  if (num < 2)
    throw VertexException("VertexKinematicConstraint::<2 states passed");
  AlgebraicMatrix jac_e(2 * num, 3);
  int num_r = 0;
  for (std::vector<KinematicState>::const_iterator i = states.begin(); i != states.end(); i++) {
    AlgebraicMatrix el_part_e(2, 3, 0);
    TrackCharge ch = i->particleCharge();
    GlobalVector mom = i->globalMomentum();
    GlobalPoint pos = i->globalPosition();
    double d_x = point.x() - pos.x();
    double d_y = point.y() - pos.y();
    double pt = mom.transverse();

    if (ch != 0) {
      //charged particle
      double a_i = -ch * i->magneticField()->inInverseGeV(pos).z();

      double pvx = mom.x() - a_i * d_y;
      double pvy = mom.y() + a_i * d_x;
      double pvt = sqrt(pvx * pvx + pvy * pvy);
      double n = a_i * (d_x * mom.x() + d_y * mom.y());
      double m = (pvx * mom.x() + pvy * mom.y());
      double k = -mom.z() / (pvt * pvt * pt * pt);

      //E jacobian matrix
      el_part_e(1, 1) = -(mom.y() + a_i * d_x);
      el_part_e(1, 2) = mom.x() - a_i * d_y;
      el_part_e(2, 1) = k * (m * mom.x() - n * mom.y());
      el_part_e(2, 2) = k * (m * mom.y() + n * mom.x());
      el_part_e(2, 3) = 1;
      jac_e.sub(2 * num_r + 1, 1, el_part_e);
    } else {
      //neutral particle
      el_part_e(1, 1) = -mom.y();
      el_part_e(1, 2) = mom.x();
      el_part_e(2, 1) = -mom.x() * mom.z() / (pt * pt);
      el_part_e(2, 2) = -mom.y() * mom.z() / (pt * pt);
      el_part_e(2, 3) = 1;
      jac_e.sub(2 * num_r + 1, 1, el_part_e);
    }
    num_r++;
  }
  return jac_e;
}

int VertexKinematicConstraint::numberOfEquations() const { return 2; }