CMSSW Reference Manual

00001 #include "RecoVertex/KinematicFit/interface/FourMomentumKinematicConstraint.h"
00002 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
00003 
00004 
00005 FourMomentumKinematicConstraint::FourMomentumKinematicConstraint(const AlgebraicVector& momentum,
00006                                                                  const AlgebraicVector& deviation)
00007 {
00008  if((momentum.num_row() != 4)||(deviation.num_row() != 4)) 
00009   throw VertexException("FourMomentumKinematicConstraint::vector of wrong size passed as 4-Momentum or Deviations");
00010  mm = momentum; 
00011  AlgebraicVector deviation_l(7,0);
00012  deviation_l(6) = momentum(3);
00013  deviation_l(5) = momentum(2);
00014  deviation_l(4) = momentum(1);
00015   
00016  double mass_sq = momentum(4)*momentum(4) - momentum(3)*momentum(3)
00017                 -momentum(2)*momentum(2) - momentum(1)*momentum(1); 
00018 
00019  if(mass_sq == 0.)throw VertexException("FourMomentumKinematicConstraint::the constraint vector passed corresponds to 0 mass");
00020 //deviation for mass calculated from deviations
00021 //of momentum
00022  deviation_l(7) = momentum(4)*momentum(4)*deviation(4)/mass_sq 
00023                 + momentum(3)*momentum(3)*deviation(3)/mass_sq
00024                 + momentum(2)*momentum(2)*deviation(2)/mass_sq
00025                 + momentum(1)*momentum(1)*deviation(1)/mass_sq;
00026 //mass sigma because of the energy 
00027  
00028  dd = deviation_l;
00029 }
00030 
00031 pair<AlgebraicVector,AlgebraicVector> FourMomentumKinematicConstraint::value(const AlgebraicVector& exPoint) const
00032 {
00033  if(exPoint.num_row() ==0 ) throw VertexException("FourMomentumKinematicConstraint::value requested for zero Linearization point");
00034 
00035 //security check for extended cartesian parametrization 
00036  int inSize = exPoint.num_row(); 
00037  if((inSize%7) !=0) throw VertexException("FourMomentumKinematicConstraint::linearization point has a wrong dimension");
00038  int nStates = inSize/7;
00039  if(nStates != 1) throw VertexException("FourMomentumKinematicConstraint::Multi particle refit is not supported in this version");
00040  AlgebraicVector pr = exPoint;
00041  AlgebraicVector vl(4,0);
00042  vl(1) += (pr(4) - mm(1));
00043  vl(2) += (pr(5) - mm(2));
00044  vl(3) += (pr(6) - mm(3));
00045  vl(7) += (sqrt(pr(4)*pr(4)+pr(5)*pr(5)+pr(6)*pr(6)+pr(7)*pr(7)) - mm(4));
00046  
00047  return pair<AlgebraicVector,AlgebraicVector>(vl,pr);
00048 }
00049  
00050 pair<AlgebraicMatrix, AlgebraicVector> FourMomentumKinematicConstraint::derivative(const AlgebraicVector& exPoint) const
00051 {
00052  if(exPoint.num_row() ==0) throw VertexException("FourMomentumKinematicConstraint::value requested for zero Linearization point");
00053 
00054 //security check for extended cartesian parametrization 
00055  int inSize = exPoint.num_row(); 
00056  if((inSize%7) !=0) throw VertexException("FourMomentumKinematicConstraint::linearization point has a wrong dimension");
00057  int nStates = inSize/7;
00058  if(nStates != 1) throw VertexException("FourMomentumKinematicConstraint::Multi particle refit is not supported in this version");
00059  AlgebraicVector pr = exPoint;
00060  AlgebraicMatrix dr(4,7,0);
00061 
00062  dr(1,4) = 1.;
00063  dr(2,5) = 1.;
00064  dr(3,6) = 1.;
00065  dr(4,7) = pr(7)/ sqrt(pr(4)*pr(4)+pr(5)*pr(5)+pr(6)*pr(6)+pr(7)*pr(7));
00066   
00067  return pair<AlgebraicMatrix,AlgebraicVector>(dr,pr);
00068 }
00069 
00070 pair<AlgebraicVector, AlgebraicVector> FourMomentumKinematicConstraint::value(const vector<RefCountedKinematicParticle> par) const
00071 {
00072  int nStates = par.size();
00073  if(nStates == 0) throw VertexException("FourMomentumKinematicConstraint::Empty vector of particles passed");
00074  if(nStates != 1) throw VertexException("FourMomentumKinematicConstraint::Multi particle refit is not supported in this version");
00075  AlgebraicVector pr = asHepVector<7>(par.front()->currentState().kinematicParameters().vector());
00076  AlgebraicVector vl(4,0);
00077  
00078  vl(1) = pr(4) - mm(1);
00079  vl(2) = pr(5) - mm(2);
00080  vl(3) = pr(6) - mm(3);
00081  vl(7) = sqrt(pr(4)*pr(4)+pr(5)*pr(5)+pr(6)*pr(6)+pr(7)*pr(7)) - mm(4);
00082  
00083  return pair<AlgebraicVector,AlgebraicVector>(vl,pr);
00084 }
00085 
00086 pair<AlgebraicMatrix, AlgebraicVector> FourMomentumKinematicConstraint::derivative(const vector<RefCountedKinematicParticle> par) const
00087 {
00088  int nStates = par.size();
00089  if(nStates == 0) throw VertexException("FourMomentumKinematicConstraint::Empty vector of particles passed");
00090  if(nStates != 1) throw VertexException("FourMomentumKinematicConstraint::Multi particle refit is not supported in this version");
00091  AlgebraicMatrix dr(4,7,0);
00092  
00093  AlgebraicVector pr = asHepVector<7>(par.front()->currentState().kinematicParameters().vector());
00094  dr(1,4) = 1.;
00095  dr(2,5) = 1.;
00096  dr(3,6) = 1.;
00097  dr(4,7) = - pr(7)/sqrt(pr(4)*pr(4)+pr(5)*pr(5)+pr(6)*pr(6)+pr(7)*pr(7));
00098  
00099  return pair<AlgebraicMatrix,AlgebraicVector>(dr,pr);
00100 }
00101 
00102 AlgebraicVector FourMomentumKinematicConstraint::deviations(int nStates) const
00103 {
00104  if(nStates == 0) throw VertexException("FourMomentumKinematicConstraint::Empty vector of particles passed");
00105  if(nStates != 1) throw VertexException("FourMomentumKinematicConstraint::Multi particle refit is not supported in this version");
00106  AlgebraicVector res = dd; 
00107  return res;
00108 }
00109 
00110 int FourMomentumKinematicConstraint::numberOfEquations() const
00111 {return 4;}
00112
FourMomentumKinematicConstraint.cc
