MultiTrackVertexLinkKinematicConstraint.cc

Go to the documentation of this file.

#include "RecoVertex/KinematicFit/interface/MultiTrackVertexLinkKinematicConstraint.h"
#include "RecoVertex/VertexPrimitives/interface/VertexException.h"

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

  //2 equations (for all tracks)
  AlgebraicVector vl(2, 0);
  double dx = point.x() - refPoint.x();
  double dy = point.y() - refPoint.y();
  double dz = point.z() - refPoint.z();
  double dT = sqrt(pow(dx, 2) + pow(dy, 2));
  double ds = sqrt(pow(dx, 2) + pow(dy, 2) + pow(dz, 2));

  double pxSum = 0, pySum = 0, pzSum = 0;
  double aSum = 0;
  for (std::vector<KinematicState>::const_iterator i = states.begin(); i != states.end(); i++) {
    double a = -i->particleCharge() * i->magneticField()->inInverseGeV(i->globalPosition()).z();
    aSum += a;

    pxSum += i->kinematicParameters()(3) - a * (point.y() - i->kinematicParameters()(1));
    pySum += i->kinematicParameters()(4) + a * (point.x() - i->kinematicParameters()(0));
    pzSum += i->kinematicParameters()(5);
  }

  double pT = sqrt(pow(pxSum, 2) + pow(pySum, 2));
  double pSum = sqrt(pow(pxSum, 2) + pow(pySum, 2) + pow(pzSum, 2));

  vl(1) = (dT - dx) / dy + (-2 * pT + sqrt(-(pow(aSum, 2) * pow(dT, 2)) + 4 * pow(pT, 2))) / (aSum * dT) +
          (-pT + pxSum) / pySum;
  vl(2) = (ds - dT) / dz + (pT - pSum) / pzSum;

  return vl;
}

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

  //2 equations (for all tracks)
  AlgebraicMatrix matrix(2, num * 7, 0);  //AlgebraicMatrix starts from 1
  double dx = point.x() - refPoint.x();
  double dy = point.y() - refPoint.y();
  double dT = sqrt(pow(dx, 2) + pow(dy, 2));

  double pxSum = 0, pySum = 0, pzSum = 0;
  double aSum = 0;
  for (std::vector<KinematicState>::const_iterator i = states.begin(); i != states.end(); i++) {
    double a = -i->particleCharge() * i->magneticField()->inInverseGeV(i->globalPosition()).z();
    aSum += a;

    pxSum += i->kinematicParameters()(3) - a * (point.y() - i->kinematicParameters()(1));
    pySum += i->kinematicParameters()(4) + a * (point.x() - i->kinematicParameters()(0));
    pzSum += i->kinematicParameters()(5);
  }

  double pT = sqrt(pow(pxSum, 2) + pow(pySum, 2));
  double pSum = sqrt(pow(pxSum, 2) + pow(pySum, 2) + pow(pzSum, 2));

  int col = 0;
  for (std::vector<KinematicState>::const_iterator i = states.begin(); i != states.end(); i++) {
    double a = -i->particleCharge() * i->magneticField()->inInverseGeV(i->globalPosition()).z();

    matrix(1, 1 + col * 7) = a * (-(pT / pow(pySum, 2)) + pxSum / pow(pySum, 2) -
                                  (4 * pySum) / (aSum * dT * sqrt(-(pow(aSum, 2) * pow(dT, 2)) + 4 * pow(pT, 2))) +
                                  (1 + (2 * pySum) / (aSum * dT)) / pT);  //dH/dx
    matrix(1, 2 + col * 7) =
        (a * (aSum * dT * (pT - pxSum) +
              2 * (-1 + (2 * pT) / sqrt(-(pow(aSum, 2) * pow(dT, 2)) + 4 * pow(pT, 2))) * pxSum * pySum)) /
        (aSum * dT * pT * pySum);  //dH/dy
    //dH/dz=0
    matrix(1, 4 + col * 7) =
        (aSum * dT * (pT - pxSum) +
         2 * (-1 + (2 * pT) / sqrt(-(pow(aSum, 2) * pow(dT, 2)) + 4 * pow(pT, 2))) * pxSum * pySum) /
        (aSum * dT * pT * pySum);  //dH/dpx
    matrix(1, 5 + col * 7) = pT / pow(pySum, 2) - pxSum / pow(pySum, 2) +
                             (4 * pySum) / (aSum * dT * sqrt(-(pow(aSum, 2) * pow(dT, 2)) + 4 * pow(pT, 2))) +
                             (-1 - (2 * pySum) / (aSum * dT)) / pT;  //dH/dpy
    //dH/dpz=0
    //dH/dm=0
    matrix(2, 1 + col * 7) = (a * (-pSum + pT) * pySum) / (pSum * pT * pzSum);  //dH/dx
    matrix(2, 2 + col * 7) = (a * (pSum - pT) * pxSum) / (pSum * pT * pzSum);   //dH/dy
    //dH/dz
    matrix(2, 4 + col * 7) = ((-(1 / pSum) + 1 / pT) * pxSum) / pzSum;   //dH/dpx
    matrix(2, 5 + col * 7) = ((-(1 / pSum) + 1 / pT) * pySum) / pzSum;   //dH/dpy
    matrix(2, 6 + col * 7) = -(1 / pSum) + (pSum - pT) / pow(pzSum, 2);  //dH/dpz
    //dH/dm=0

    col++;
  }

  return matrix;
}

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

  //2 equations (for all tracks)
  AlgebraicMatrix matrix(2, 3, 0);
  double dx = point.x() - refPoint.x();
  double dy = point.y() - refPoint.y();
  double dz = point.z() - refPoint.z();
  double dT = sqrt(pow(dx, 2) + pow(dy, 2));
  double ds = sqrt(pow(dx, 2) + pow(dy, 2) + pow(dz, 2));

  double pxSum = 0, pySum = 0, pzSum = 0, aSum = 0;
  for (std::vector<KinematicState>::const_iterator i = states.begin(); i != states.end(); i++) {
    double a = -i->particleCharge() * i->magneticField()->inInverseGeV(i->globalPosition()).z();
    aSum += a;

    pxSum += i->kinematicParameters()(3) - a * (point.y() - i->kinematicParameters()(1));
    pySum += i->kinematicParameters()(4) + a * (point.x() - i->kinematicParameters()(0));
    pzSum += i->kinematicParameters()(5);
  }
  double pT = sqrt(pow(pxSum, 2) + pow(pySum, 2));
  double pSum = sqrt(pow(pxSum, 2) + pow(pySum, 2) + pow(pzSum, 2));

  matrix(1, 1) =
      (-1 + dx / dT) / dy +
      (2 * dx * pT * (1 - (2 * pT) / sqrt(-(pow(aSum, 2) * pow(dT, 2)) + 4 * pow(pT, 2)))) / (aSum * pow(dT, 3)) +
      aSum * (-(1 / pT) + pT / pow(pySum, 2) - pxSum / pow(pySum, 2)) +
      (2 * (-(1 / pT) + 2 / sqrt(-(pow(aSum, 2) * pow(dT, 2)) + 4 * pow(pT, 2))) * pySum) / dT;  //dH/dxv
  matrix(1, 2) = 1 / dT + (-dT + dx) / pow(dy, 2) -
                 (dy * (-2 * pT + sqrt(-(pow(aSum, 2) * pow(dT, 2)) + 4 * pow(pT, 2)))) / (aSum * pow(dT, 3)) -
                 ((-2 + sqrt(4 - (pow(aSum, 2) * pow(dT, 2)) / pow(pT, 2))) * pxSum) / (dT * pT) -
                 (aSum * (dy * pow(pT, 2) + aSum * pow(dT, 2) * pxSum)) /
                     (dT * pow(pT, 2) * sqrt(-(pow(aSum, 2) * pow(dT, 2)) + 4 * pow(pT, 2))) +
                 (aSum * (-pT + pxSum)) / (pT * pySum);  //dH/dyv
  //dH/dzv=0
  matrix(2, 1) = ((1 / ds - 1 / dT) * dx) / dz + (aSum * (pSum - pT) * pySum) / (pSum * pT * pzSum);  //dH/dxv
  matrix(2, 2) = ((1 / ds - 1 / dT) * dy) / dz - (aSum * (pSum - pT) * pxSum) / (pSum * pT * pzSum);  //dH/dyv
  matrix(2, 3) = 1 / ds + (-ds + dT) / pow(dz, 2);                                                    //dH/dzv

  return matrix;
}

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