CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
FourMomentumKinematicConstraint.cc
Go to the documentation of this file.
3 
4 
6  const AlgebraicVector& deviation)
7 {
8  if((momentum.num_row() != 4)||(deviation.num_row() != 4))
9  throw VertexException("FourMomentumKinematicConstraint::vector of wrong size passed as 4-Momentum or Deviations");
10  mm = momentum;
11  AlgebraicVector deviation_l(7,0);
12  deviation_l(6) = momentum(3);
13  deviation_l(5) = momentum(2);
14  deviation_l(4) = momentum(1);
15 
16  double mass_sq = momentum(4)*momentum(4) - momentum(3)*momentum(3)
17  -momentum(2)*momentum(2) - momentum(1)*momentum(1);
18 
19  if(mass_sq == 0.)throw VertexException("FourMomentumKinematicConstraint::the constraint vector passed corresponds to 0 mass");
20 //deviation for mass calculated from deviations
21 //of momentum
22  deviation_l(7) = momentum(4)*momentum(4)*deviation(4)/mass_sq
23  + momentum(3)*momentum(3)*deviation(3)/mass_sq
24  + momentum(2)*momentum(2)*deviation(2)/mass_sq
25  + momentum(1)*momentum(1)*deviation(1)/mass_sq;
26 //mass sigma because of the energy
27 
28  dd = deviation_l;
29 }
30 
31 std::pair<AlgebraicVector,AlgebraicVector> FourMomentumKinematicConstraint::value(const AlgebraicVector& exPoint) const
32 {
33  if(exPoint.num_row() ==0 ) throw VertexException("FourMomentumKinematicConstraint::value requested for zero Linearization point");
34 
35 //security check for extended cartesian parametrization
36  int inSize = exPoint.num_row();
37  if((inSize%7) !=0) throw VertexException("FourMomentumKinematicConstraint::linearization point has a wrong dimension");
38  int nStates = inSize/7;
39  if(nStates != 1) throw VertexException("FourMomentumKinematicConstraint::Multi particle refit is not supported in this version");
40  AlgebraicVector pr = exPoint;
41  AlgebraicVector vl(4,0);
42  vl(1) += (pr(4) - mm(1));
43  vl(2) += (pr(5) - mm(2));
44  vl(3) += (pr(6) - mm(3));
45  vl(7) += (sqrt(pr(4)*pr(4)+pr(5)*pr(5)+pr(6)*pr(6)+pr(7)*pr(7)) - mm(4));
46 
47  return std::pair<AlgebraicVector,AlgebraicVector>(vl,pr);
48 }
49 
50 std::pair<AlgebraicMatrix, AlgebraicVector> FourMomentumKinematicConstraint::derivative(const AlgebraicVector& exPoint) const
51 {
52  if(exPoint.num_row() ==0) throw VertexException("FourMomentumKinematicConstraint::value requested for zero Linearization point");
53 
54 //security check for extended cartesian parametrization
55  int inSize = exPoint.num_row();
56  if((inSize%7) !=0) throw VertexException("FourMomentumKinematicConstraint::linearization point has a wrong dimension");
57  int nStates = inSize/7;
58  if(nStates != 1) throw VertexException("FourMomentumKinematicConstraint::Multi particle refit is not supported in this version");
59  AlgebraicVector pr = exPoint;
60  AlgebraicMatrix dr(4,7,0);
61 
62  dr(1,4) = 1.;
63  dr(2,5) = 1.;
64  dr(3,6) = 1.;
65  dr(4,7) = pr(7)/ sqrt(pr(4)*pr(4)+pr(5)*pr(5)+pr(6)*pr(6)+pr(7)*pr(7));
66 
67  return std::pair<AlgebraicMatrix,AlgebraicVector>(dr,pr);
68 }
69 
70 std::pair<AlgebraicVector, AlgebraicVector> FourMomentumKinematicConstraint::value(const std::vector<RefCountedKinematicParticle> par) const
71 {
72  int nStates = par.size();
73  if(nStates == 0) throw VertexException("FourMomentumKinematicConstraint::Empty vector of particles passed");
74  if(nStates != 1) throw VertexException("FourMomentumKinematicConstraint::Multi particle refit is not supported in this version");
75  AlgebraicVector pr = asHepVector<7>(par.front()->currentState().kinematicParameters().vector());
76  AlgebraicVector vl(4,0);
77 
78  vl(1) = pr(4) - mm(1);
79  vl(2) = pr(5) - mm(2);
80  vl(3) = pr(6) - mm(3);
81  vl(7) = sqrt(pr(4)*pr(4)+pr(5)*pr(5)+pr(6)*pr(6)+pr(7)*pr(7)) - mm(4);
82 
83  return std::pair<AlgebraicVector,AlgebraicVector>(vl,pr);
84 }
85 
86 std::pair<AlgebraicMatrix, AlgebraicVector> FourMomentumKinematicConstraint::derivative(const std::vector<RefCountedKinematicParticle> par) const
87 {
88  int nStates = par.size();
89  if(nStates == 0) throw VertexException("FourMomentumKinematicConstraint::Empty vector of particles passed");
90  if(nStates != 1) throw VertexException("FourMomentumKinematicConstraint::Multi particle refit is not supported in this version");
91  AlgebraicMatrix dr(4,7,0);
92 
93  AlgebraicVector pr = asHepVector<7>(par.front()->currentState().kinematicParameters().vector());
94  dr(1,4) = 1.;
95  dr(2,5) = 1.;
96  dr(3,6) = 1.;
97  dr(4,7) = - pr(7)/sqrt(pr(4)*pr(4)+pr(5)*pr(5)+pr(6)*pr(6)+pr(7)*pr(7));
98 
99  return std::pair<AlgebraicMatrix,AlgebraicVector>(dr,pr);
100 }
101 
103 {
104  if(nStates == 0) throw VertexException("FourMomentumKinematicConstraint::Empty vector of particles passed");
105  if(nStates != 1) throw VertexException("FourMomentumKinematicConstraint::Multi particle refit is not supported in this version");
106  AlgebraicVector res = dd;
107  return res;
108 }
109 
111 {return 4;}
112 
virtual std::pair< AlgebraicMatrix, AlgebraicVector > derivative(const AlgebraicVector &exPoint) const
Common base class.
virtual std::pair< AlgebraicVector, AlgebraicVector > value(const AlgebraicVector &exPoint) const
CLHEP::HepMatrix AlgebraicMatrix
T sqrt(T t)
Definition: SSEVec.h:46
virtual AlgebraicVector deviations(int nStates) const
CLHEP::HepVector AlgebraicVector
FourMomentumKinematicConstraint(const AlgebraicVector &momentum, const AlgebraicVector &deviation)