FourMomentumKinematicConstraint.cc

Go to the documentation of this file.

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


FourMomentumKinematicConstraint::FourMomentumKinematicConstraint(const AlgebraicVector& momentum,
                                                                 const AlgebraicVector& deviation)
{
 if((momentum.num_row() != 4)||(deviation.num_row() != 4)) 
  throw VertexException("FourMomentumKinematicConstraint::vector of wrong size passed as 4-Momentum or Deviations");
 mm = momentum; 
 AlgebraicVector deviation_l(7,0);
 deviation_l(6) = momentum(3);
 deviation_l(5) = momentum(2);
 deviation_l(4) = momentum(1);
  
 double mass_sq = momentum(4)*momentum(4) - momentum(3)*momentum(3)
                -momentum(2)*momentum(2) - momentum(1)*momentum(1); 

 if(mass_sq == 0.)throw VertexException("FourMomentumKinematicConstraint::the constraint vector passed corresponds to 0 mass");
//deviation for mass calculated from deviations
//of momentum
 deviation_l(7) = momentum(4)*momentum(4)*deviation(4)/mass_sq 
                + momentum(3)*momentum(3)*deviation(3)/mass_sq
        + momentum(2)*momentum(2)*deviation(2)/mass_sq
        + momentum(1)*momentum(1)*deviation(1)/mass_sq;
//mass sigma because of the energy 
 
 dd = deviation_l;
}

std::pair<AlgebraicVector,AlgebraicVector> FourMomentumKinematicConstraint::value(const AlgebraicVector& exPoint) const
{
 if(exPoint.num_row() ==0 ) throw VertexException("FourMomentumKinematicConstraint::value requested for zero Linearization point");

//security check for extended cartesian parametrization 
 int inSize = exPoint.num_row(); 
 if((inSize%7) !=0) throw VertexException("FourMomentumKinematicConstraint::linearization point has a wrong dimension");
 int nStates = inSize/7;
 if(nStates != 1) throw VertexException("FourMomentumKinematicConstraint::Multi particle refit is not supported in this version");
 AlgebraicVector pr = exPoint;
 AlgebraicVector vl(4,0);
 vl(1) += (pr(4) - mm(1));
 vl(2) += (pr(5) - mm(2));
 vl(3) += (pr(6) - mm(3));
 vl(7) += (sqrt(pr(4)*pr(4)+pr(5)*pr(5)+pr(6)*pr(6)+pr(7)*pr(7)) - mm(4));
 
 return std::pair<AlgebraicVector,AlgebraicVector>(vl,pr);
}
 
std::pair<AlgebraicMatrix, AlgebraicVector> FourMomentumKinematicConstraint::derivative(const AlgebraicVector& exPoint) const
{
 if(exPoint.num_row() ==0) throw VertexException("FourMomentumKinematicConstraint::value requested for zero Linearization point");

//security check for extended cartesian parametrization 
 int inSize = exPoint.num_row(); 
 if((inSize%7) !=0) throw VertexException("FourMomentumKinematicConstraint::linearization point has a wrong dimension");
 int nStates = inSize/7;
 if(nStates != 1) throw VertexException("FourMomentumKinematicConstraint::Multi particle refit is not supported in this version");
 AlgebraicVector pr = exPoint;
 AlgebraicMatrix dr(4,7,0);

 dr(1,4) = 1.;
 dr(2,5) = 1.;
 dr(3,6) = 1.;
 dr(4,7) = pr(7)/ sqrt(pr(4)*pr(4)+pr(5)*pr(5)+pr(6)*pr(6)+pr(7)*pr(7));
  
 return std::pair<AlgebraicMatrix,AlgebraicVector>(dr,pr);
}

std::pair<AlgebraicVector, AlgebraicVector> FourMomentumKinematicConstraint::value(const std::vector<RefCountedKinematicParticle> &par) const
{
 int nStates = par.size();
 if(nStates == 0) throw VertexException("FourMomentumKinematicConstraint::Empty vector of particles passed");
 if(nStates != 1) throw VertexException("FourMomentumKinematicConstraint::Multi particle refit is not supported in this version");
 AlgebraicVector pr = asHepVector<7>(par.front()->currentState().kinematicParameters().vector());
 AlgebraicVector vl(4,0);
 
 vl(1) = pr(4) - mm(1);
 vl(2) = pr(5) - mm(2);
 vl(3) = pr(6) - mm(3);
 vl(7) = sqrt(pr(4)*pr(4)+pr(5)*pr(5)+pr(6)*pr(6)+pr(7)*pr(7)) - mm(4);
 
 return std::pair<AlgebraicVector,AlgebraicVector>(vl,pr);
}

std::pair<AlgebraicMatrix, AlgebraicVector> FourMomentumKinematicConstraint::derivative(const std::vector<RefCountedKinematicParticle> &par) const
{
 int nStates = par.size();
 if(nStates == 0) throw VertexException("FourMomentumKinematicConstraint::Empty vector of particles passed");
 if(nStates != 1) throw VertexException("FourMomentumKinematicConstraint::Multi particle refit is not supported in this version");
 AlgebraicMatrix dr(4,7,0);
 
 AlgebraicVector pr = asHepVector<7>(par.front()->currentState().kinematicParameters().vector());
 dr(1,4) = 1.;
 dr(2,5) = 1.;
 dr(3,6) = 1.;
 dr(4,7) = - pr(7)/sqrt(pr(4)*pr(4)+pr(5)*pr(5)+pr(6)*pr(6)+pr(7)*pr(7));
 
 return std::pair<AlgebraicMatrix,AlgebraicVector>(dr,pr);
}

AlgebraicVector FourMomentumKinematicConstraint::deviations(int nStates) const
{
 if(nStates == 0) throw VertexException("FourMomentumKinematicConstraint::Empty vector of particles passed");
 if(nStates != 1) throw VertexException("FourMomentumKinematicConstraint::Multi particle refit is not supported in this version");
 AlgebraicVector res = dd; 
 return res;
}

int FourMomentumKinematicConstraint::numberOfEquations() const
{return 4;}