VertexKinematicConstraintT.cc

Go to the documentation of this file.

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

VertexKinematicConstraintT::VertexKinematicConstraintT()
{}

VertexKinematicConstraintT::~VertexKinematicConstraintT()
{}


void VertexKinematicConstraintT::init(const std::vector<KinematicState>& states,
                      const GlobalPoint& ipoint,  const GlobalVector& fieldValue) {
  int num = states.size();
  if(num!=2) throw VertexException("VertexKinematicConstraintT !=2 states passed");
  double mfz = fieldValue.z();
  
  int j=0;
  for(std::vector<KinematicState>::const_iterator i = states.begin(); i != states.end(); i++)
    {
      mom[j] = i->globalMomentum();
      dpos[j] = ipoint - i->globalPosition();
      a_i[j] = - i->particleCharge() * mfz;

      double pvx = mom[j].x() - a_i[j]*dpos[j].y();
      double pvy = mom[j].y() + a_i[j]*dpos[j].x();
      double pvt2 = pvx*pvx+pvy*pvy;
      novera[j] = (dpos[j].x() * mom[j].x() + dpos[j].y()*mom[j].y());
      n[j] = a_i[j]*novera[j];
      m[j] = (pvx*mom[j].x() + pvy*mom[j].y());
      k[j] = -mom[j].z()/(mom[j].perp2()*pvt2);
      delta[j] = std::atan2(n[j],m[j]);

      ++j;
    }
}



void VertexKinematicConstraintT::fillValue() const
{
  //it is 2 equations per track
  for(int j = 0; j!=2; ++j) {

    if(a_i[j] !=0) {
      //vector of values
      super::vl(j*2) = dpos[j].y()*mom[j].x() - dpos[j].x()*mom[j].y() -a_i[j]*(dpos[j].x()*dpos[j].x() + dpos[j].y()*dpos[j].y())*0.5;
      super::vl(j*2 +1) = dpos[j].z() - mom[j].z()*delta[j]/a_i[j];
    }else{      
      //neutral particle
      double pt2Inverse = 1./mom[j].perp2();
      super::vl(j*2) = dpos[j].y()*mom[j].x() - dpos[j].x()*mom[j].y();
      super::vl(j*2 +1) = dpos[j].z() - mom[j].z()*(dpos[j].x() * mom[j].x() + dpos[j].y() * mom[j].y())*pt2Inverse;
    }
  }
}

void VertexKinematicConstraintT::fillParametersDerivative() const {
  ROOT::Math::SMatrix<double,2,7> el_part_d;
  for(int j = 0; j!=2; ++j) {
    
    if(a_i[j] !=0) {
      
      //charged particle
            
      //D Jacobian matrix
      el_part_d(0,0) =  mom[j].y() + a_i[j]*dpos[j].x();
      el_part_d(0,1) = -mom[j].x() + a_i[j]*dpos[j].y();
      el_part_d(1,0) =  -k[j]*(m[j]*mom[j].x() - n[j]*mom[j].y());
      el_part_d(1,1) =  -k[j]*(m[j]*mom[j].y() + n[j]*mom[j].x());
      el_part_d(1,2) = -1.;
      el_part_d(0,3) = dpos[j].y();
      el_part_d(0,4) = -dpos[j].x();
      el_part_d(1,3) = k[j]*(m[j]*dpos[j].x() - novera[j]*(2*mom[j].x() - a_i[j]*dpos[j].y()));
      el_part_d(1,4) = k[j]*(m[j]*dpos[j].y() - novera[j]*(2*mom[j].y() + a_i[j]*dpos[j].x()));
      el_part_d(1,5) = -delta[j] /a_i[j];
      super::jac_d().Place_at(el_part_d,j*2, j*7);
    }else{
      //neutral particle
      double pt2Inverse = 1./mom[j].perp2();
      el_part_d(0,0) =  mom[j].y();
      el_part_d(0,1) = -mom[j].x();
      el_part_d(1,0) =  mom[j].x() * (mom[j].z()*pt2Inverse);
      el_part_d(1,1) =  mom[j].y() * (mom[j].z()*pt2Inverse);
      el_part_d(1,2) = -1.;
      el_part_d(0,3) = dpos[j].y();
      el_part_d(0,4) = -dpos[j].x();
      el_part_d(1,3) = 2*(dpos[j].x()*mom[j].x()+dpos[j].y()*mom[j].y())*pt2Inverse*mom[j].x()*(mom[j].z()*pt2Inverse) - dpos[j].x()*(mom[j].z()*pt2Inverse);
      el_part_d(1,4) = 2*(dpos[j].x()*mom[j].x()+dpos[j].y()*mom[j].y())*pt2Inverse*mom[j].y()*(mom[j].z()*pt2Inverse) - dpos[j].x()*(mom[j].z()*pt2Inverse);
      el_part_d(1,5) = - (dpos[j].x()*mom[j].x()+dpos[j].y()*mom[j].y())*pt2Inverse;
      super::jac_d().Place_at(el_part_d,j*2, j*7);
    }
  }
}

void VertexKinematicConstraintT::fillPositionDerivative() const {

  ROOT::Math::SMatrix<double,2,3>  el_part_e;
  for(int j = 0; j!=2; ++j) {

    if(a_i[j] !=0) {
      
      //charged particle

      //E jacobian matrix
      el_part_e(0,0) = -(mom[j].y() + a_i[j]*dpos[j].x());
      el_part_e(0,1) =   mom[j].x() - a_i[j]*dpos[j].y();
      el_part_e(1,0) = k[j]*(m[j]*mom[j].x() - n[j]*mom[j].y());
      el_part_e(1,1) = k[j]*(m[j]*mom[j].y() + n[j]*mom[j].x());
      el_part_e(1,2) = 1;
      super::jac_e().Place_at(el_part_e,2*j,0);
    }else{      
      //neutral particle
      double pt2Inverse = 1./mom[j].perp2();
      el_part_e(0,0) = - mom[j].y();
      el_part_e(0,1) = mom[j].x();
      el_part_e(1,0) = -mom[j].x()*mom[j].z()*pt2Inverse;
      el_part_e(1,1) = -mom[j].y()*mom[j].z()*pt2Inverse;
      el_part_e(1,2) = 1;
      super::jac_e().Place_at(el_part_e,2*j,0);
    }
  }
}

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