/data/refman/pasoursint/CMSSW_4_1_8_patch13/src/RecoVertex/KinematicFit/src/MomentumKinematicConstraint.cc

Go to the documentation of this file.

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


MomentumKinematicConstraint::MomentumKinematicConstraint(const AlgebraicVector& momentum,
                                                         const AlgebraicVector& dev)
{
 if((momentum.num_row() != 3) || (dev.num_row() != 3))
     throw VertexException("MomentumKinemticConstraint::Momentum or Deviation vector passed is not 3-dimensional");
 mm = momentum;
 AlgebraicVector dev_l(7,0);
 dev_l(4) = dev(1) * dev(1);
 dev_l(5) = dev(2) * dev(2);
 dev_l(6) = dev(3) * dev(3);
 dd = dev_l;
}

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

//security check for extended cartesian parametrization 
 int inSize = exPoint.num_row(); 
 if((inSize%7) !=0) throw VertexException("MomentumKinematicConstraint::linearization point has a wrong dimension");
 int nStates = inSize/7;
 if(nStates != 1) throw VertexException("MomentumKinematicConstraint::Multistate refit is not foreseen for this constraint"); 
 AlgebraicVector pr = exPoint;
 AlgebraicVector vl(3,0);
 vl(1) = pr(4) - mm(1);
 vl(2) = pr(5) - mm(2);
 vl(3) = pr(6) - mm(3);
 return std::pair<AlgebraicVector,AlgebraicVector>(vl,pr); 
}

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

//security check for extended cartesian parametrization 
 int inSize = exPoint.num_row(); 
 if((inSize%7) !=0) throw VertexException("MomentumKinematicConstraint::linearization point has a wrong dimension");
 int nStates = inSize/7;
 if(nStates != 1) throw VertexException("MomentumKinematicConstraint::Multistate refit is not foreseen for this constraint"); 

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

std::pair<AlgebraicVector, AlgebraicVector> MomentumKinematicConstraint::value(const std::vector<RefCountedKinematicParticle> par) const
{
 int nStates = par.size(); 
 if(nStates == 0) throw VertexException("MomentumKinematicConstraint::Empty vector of particles passed");
 if(nStates != 1) throw VertexException("MomentumKinematicConstraint::Multistate refit is not foreseen for this constraint");  
 AlgebraicVector point = asHepVector<7>(par.front()->currentState().kinematicParameters().vector());
 AlgebraicVector vl(3,0);
 vl(1) = point(4) - mm(1);
 vl(2) = point(5) - mm(2);
 vl(3) = point(6) - mm(3);
 return std::pair<AlgebraicVector,AlgebraicVector>(vl,point);
}

std::pair<AlgebraicMatrix, AlgebraicVector> MomentumKinematicConstraint::derivative(const std::vector<RefCountedKinematicParticle> par) const
{
 int nStates = par.size();
 if(nStates == 0) throw VertexException("MomentumKinematicConstraint::Empty vector of particles passed");
 if(nStates != 1) throw VertexException("MomentumKinematicConstraint::Multistate refit is not foreseen for this constraint");
 AlgebraicVector point = asHepVector<7>(par.front()->currentState().kinematicParameters().vector());
 AlgebraicMatrix dr(3,7,0);
 dr(1,4) = 1.;
 dr(2,5) = 1.;
 dr(3,6) = 1.;
 return std::pair<AlgebraicMatrix,AlgebraicVector>(dr,point);
}

AlgebraicVector MomentumKinematicConstraint::deviations(int nStates) const
{
 if(nStates == 0) throw VertexException("MomentumKinematicConstraint::Empty vector of particles passed");
 if(nStates != 1) throw VertexException("MomentumKinematicConstraint::Multistate refit is not foreseen for this constraint");
 AlgebraicVector res = dd;
 return res;
}

int MomentumKinematicConstraint::numberOfEquations() const
{return 3;}