ColinearityKinematicConstraint.cc

Go to the documentation of this file.

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

ColinearityKinematicConstraint::ColinearityKinematicConstraint(ConstraintDim dim)
{
  dimension = dim;
  if (dimension == Phi) size = 1;
  else size = 2;
}

AlgebraicVector  ColinearityKinematicConstraint::value(const std::vector<KinematicState> states,
                        const GlobalPoint& point) const
{
  if(states.size()<2) throw VertexException("ColinearityKinematicConstraint::<2 states passed");
  AlgebraicVector res(size,0);

  double a_1 = -states[0].particleCharge()*states[0].magneticField()->inInverseGeV(states[0].globalPosition()).z();
  double a_2 = -states[1].particleCharge()*states[1].magneticField()->inInverseGeV(states[1].globalPosition()).z();
  
  AlgebraicVector7 p1 = states[0].kinematicParameters().vector();
  AlgebraicVector7 p2 = states[1].kinematicParameters().vector();

  double p1vx = p1(3) - a_1*(point.y() - p1(1));
  double p1vy = p1(4) + a_1*(point.x() - p1(0));
  double pt1  = sqrt(p1(3)*p1(3)+p1(4)*p1(4));

  double p2vx = p2(3) - a_2*(point.y() - p2(1));
  double p2vy = p2(4) + a_2*(point.x() - p2(0));
  double pt2  = sqrt(p2(3)*p2(3)+p2(4)*p2(4));

  // H_phi:
  res(1)  = atan2(p1vy,p1vx) - atan2(p2vy,p2vx);
  if ( res(1) >  M_PI ) res(1) -= 2.0*M_PI;
  if ( res(1) <= -M_PI ) res(1) += 2.0*M_PI;
  // H_theta:
  if (dimension == PhiTheta) {  
    res(2)  = atan2(pt1,p1(5)) - atan2(pt2,p2(5));
    if ( res(2) >  M_PI ) res(2) -= 2.0*M_PI;
    if ( res(2) <= -M_PI ) res(2) += 2.0*M_PI;
  }

// cout << res(1) << " "<<res(2) << "\n ";// res(2)

  return res;
}

AlgebraicMatrix ColinearityKinematicConstraint::parametersDerivative(const std::vector<KinematicState> states,
                                      const GlobalPoint& point) const
{
  int n_st = states.size();
  if(n_st<2) throw VertexException("ColinearityKinematicConstraint::<2 states passed");
  AlgebraicMatrix res(size,n_st*7,0);

  double a_1 = -states[0].particleCharge()*states[0].magneticField()->inInverseGeV(states[0].globalPosition()).z();
  double a_2 = -states[1].particleCharge()*states[1].magneticField()->inInverseGeV(states[1].globalPosition()).z();

  AlgebraicVector7 p1 = states[0].kinematicParameters().vector();
  AlgebraicVector7 p2 = states[1].kinematicParameters().vector();

  double p1vx = p1(3) - a_1*(point.y() - p1(1));
  double p1vy = p1(4) + a_1*(point.x() - p1(0));
  double k1 = 1.0/(p1vx*p1vx + p1vy*p1vy);
  double pt1 = sqrt(p1(3)*p1(3)+p1(4)*p1(4));
  double pTot1  = sqrt(p1(3)*p1(3)+p1(4)*p1(4)+p1(5)*p1(5));

  double p2vx = p2(3) - a_2*(point.y() - p2(1));
  double p2vy = p2(4) + a_2*(point.x() - p2(0));
  double k2 = 1.0/(p2vx*p2vx + p2vy*p2vy);
  double pt2  = sqrt(p2(3)*p2(3)+p2(4)*p2(4));
  double pTot2   = sqrt(p2(3)*p2(3)+p2(4)*p2(4)+p2(5)*p2(5));

  // H_phi:

  //x1 and x2 derivatives: 1st and 8th elements
  res(1,1) =  -k1*a_1*p1vx;
  res(1,8) =   k2*a_2*p2vx;

  //y1 and y2 derivatives: 2nd and 9th elements:
  res(1,2) = -k1*a_1*p1vy;
  res(1,9) =  k2*a_2*p2vy;

  //z1 and z2 components: 3d and 10th elmnets stay 0:
  res(1,3)  = 0.; res(1,10) = 0.;

  //px1 and px2 components: 4th and 11th elements:
  res(1,4)  = -k1*p1vy;
  res(1,11) =  k2*p2vy;

  //py1 and py2 components: 5th and 12 elements:
  res(1,5)  =  k1*p1vx;
  res(1,12) = -k2*p2vx;


  //pz1 and pz2 components: 6th and 13 elements:
  res(1,6)  = 0.; res(1,13) = 0.;
  //mass components: 7th and 14th elements:
  res(1,7)  = 0.; res(1,14) = 0.;

  if (dimension == PhiTheta)  {
    // H_theta:
    //x1 and x2 derivatives: 1st and 8th elements
    res(2,1) =  0.; res(2,8) = 0.;

    //y1 and y2 derivatives: 2nd and 9th elements:
    res(2,2) = 0.; res(2,9) = 0.;

    //z1 and z2 components: 3d and 10th elmnets stay 0:
    res(2,3) = 0.; res(2,10) = 0.;

    res(2,4)  =   p1(5)*p1(3) / (pTot1*pTot1*pt1);
    res(2,11) = - p2(5)*p2(3) / (pTot2*pTot2*pt2);

    //py1 and py2 components: 5th and 12 elements:
    res(2,5)  =   p1(5)*p1(4) / (pTot1*pTot1*pt1);
    res(2,12) = - p2(5)*p2(4) / (pTot2*pTot2*pt2);

    //pz1 and pz2 components: 6th and 13 elements:
    res(2,6)  = - pt1/(pTot1*pTot1);
    res(2,13) =   pt2/(pTot2*pTot2);
    //mass components: 7th and 14th elements:
    res(2,7)  = 0.; res(2,14) = 0.;
  }

  return res;
}

AlgebraicMatrix ColinearityKinematicConstraint::positionDerivative(const std::vector<KinematicState> states,
                                    const GlobalPoint& point) const
{
  AlgebraicMatrix res(size,3,0);
  if(states.size()<2) throw VertexException("ColinearityKinematicConstraint::<2 states passed");

  double a_1 = -states[0].particleCharge() * states[0].magneticField()->inInverseGeV(states[0].globalPosition()).z();
  double a_2 = -states[1].particleCharge() * states[1].magneticField()->inInverseGeV(states[1].globalPosition()).z();

  AlgebraicVector7 p1 = states[0].kinematicParameters().vector();
  AlgebraicVector7 p2 = states[1].kinematicParameters().vector();

  double p1vx = p1(3) - a_1*(point.y() - p1(1));
  double p1vy = p1(4) + a_1*(point.x() - p1(0));
  double k1 = 1.0/(p1vx*p1vx + p1vy*p1vy);
  //double pt1 = sqrt(p1(3)*p1(3)+p1(4)*p1(4));

  double p2vx = p2(3) - a_2*(point.y() - p2(1));
  double p2vy = p2(4) + a_2*(point.x() - p2(0));
  double k2 = 1.0/(p2vx*p2vx + p2vy*p2vy);
  //double pt2  = sqrt(p2(3)*p2(3)+p2(4)*p2(4));

  // H_phi:

  // xv component
  res(1,1) = k1*a_1*p1vx - k2*a_2*p2vx;

  //yv component
  res(1,2) = k1*a_1*p1vy - k2*a_2*p2vy;

  //zv component
  res(1,3) = 0.;

  // H_theta: no correlation with vertex position
  if (dimension == PhiTheta) {
    res(2,1) = 0.;
    res(2,2) = 0.;
    res(2,3) = 0.;
  }

  return res;
}