CMS 3D CMS Logo
SmartPointingConstraint.cc
Go to the documentation of this file.
1 #include "RecoVertex/KinematicFit/interface/SmartPointingConstraint.h"
2 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
3 
4 std::pair<AlgebraicVector, AlgebraicVector> SmartPointingConstraint::value(const AlgebraicVector& exPoint) const {
5  if (exPoint.num_row() == 0)
6  throw VertexException("PointingKinematicConstraint::value requested for zero Linearization point");
7 
8  //security check for extended cartesian parametrization
9  int inSize = exPoint.num_row();
10  if ((inSize % 7) != 0)
11  throw VertexException("PointingKinematicConstraint::linearization point has a wrong dimension");
12  int nStates = inSize / 7;
13  if (nStates != 1)
14  throw VertexException("PointingKinematicConstraint::Current version does not support the multistate refit");
15 
16  AlgebraicVector lPar = exPoint;
17  AlgebraicVector vl(2, 0);
18 
19  //vector of values 1x2 for given particle
20  AlgebraicVector lValue = makeValue(lPar).first;
21  vl(1) = lValue(1);
22  vl(2) = lValue(2);
23  return std::pair<AlgebraicVector, AlgebraicVector>(vl, lPar);
24 }
25 
26 std::pair<AlgebraicMatrix, AlgebraicVector> SmartPointingConstraint::derivative(const AlgebraicVector& exPoint) const {
27  if (exPoint.num_row() == 0)
28  throw VertexException("PointingKinematicConstraint::value requested for zero Linearization point");
29 
30  //security check for extended cartesian parametrization
31  int inSize = exPoint.num_row();
32  if ((inSize % 7) != 0)
33  throw VertexException("PointingKinematicConstraint::linearization point has a wrong dimension");
34  int nStates = inSize / 7;
35  if (nStates != 1)
36  throw VertexException("PointingKinematicConstraint::Current version does not support the multistate refit");
37  AlgebraicVector lPar = exPoint;
38 
39  //2x7 derivative matrix for given particle
40  AlgebraicMatrix lDeriv = makeDerivative(lPar).first;
41  AlgebraicMatrix dr(2, 7, 0);
42  dr.sub(1, 1, lDeriv);
43  return std::pair<AlgebraicMatrix, AlgebraicVector>(dr, lPar);
44 }
45 
46 std::pair<AlgebraicMatrix, AlgebraicVector> SmartPointingConstraint::derivative(
47  const std::vector<RefCountedKinematicParticle>& par) const {
48  int nStates = par.size();
49  if (nStates == 0)
50  throw VertexException("PointingKinematicConstraint::Empty vector of particles passed");
51  if (nStates != 1)
52  throw VertexException("PointingKinematicConstraint::Current version does not support the multistate refit");
53 
54  AlgebraicMatrix dr(2, 7, 0);
55  AlgebraicVector lPoint = asHepVector<7>(par.front()->currentState().kinematicParameters().vector());
56 
57  //2x7 derivative matrix for given state
58  AlgebraicMatrix lDeriv = makeDerivative(lPoint).first;
59  dr.sub(1, 1, lDeriv);
60  // cout<<"Derivative returned: "<<dr<<endl;
61  // cout<<"For the value: "<<lPoint<<endl;
62  return std::pair<AlgebraicMatrix, AlgebraicVector>(dr, lPoint);
63 }
64 
65 std::pair<AlgebraicVector, AlgebraicVector> SmartPointingConstraint::value(
66  const std::vector<RefCountedKinematicParticle>& par) const {
67  int nStates = par.size();
68  if (nStates == 0)
69  throw VertexException("PointingKinematicConstraint::Empty vector of particles passed");
70  if (nStates != 1)
71  throw VertexException("PointingKinematicConstraint::Current version does not support the multistate refit");
72  AlgebraicVector vl(2, 0);
73  AlgebraicVector lPoint = asHepVector<7>(par.front()->currentState().kinematicParameters().vector());
74  vl(1) = makeValue(lPoint).first(1);
75  vl(2) = makeValue(lPoint).first(2);
76  // cout<<"Value returned: "<<vl<<endl;
77  // cout<<"For the point: "<<lPoint<<endl;
78 
79  return std::pair<AlgebraicVector, AlgebraicVector>(vl, lPoint);
80 }
81 
82 AlgebraicVector SmartPointingConstraint::deviations(int nStates) const { return AlgebraicVector(7 * nStates, 0); }
83 
84 int SmartPointingConstraint::numberOfEquations() const { return 2; }
85 
86 std::pair<AlgebraicVector, AlgebraicVector> SmartPointingConstraint::makeValue(const AlgebraicVector& exPoint) const {
87  // cout<<"Make value called"<<endl;
88  AlgebraicVector vl(2, 0);
89  AlgebraicVector point = exPoint;
90  double dx = point(1) - refPoint.x();
91  double dy = point(2) - refPoint.y();
92  double dz = point(3) - refPoint.z();
93  double px = point(4);
94  double py = point(5);
95  double pz = point(6);
96 
97  //full angle solution: sin(alpha - betha) = 0
98  //sign swap allowed
99  double cos_phi_p = px / sqrt(px * px + py * py);
100  double sin_phi_p = py / sqrt(px * px + py * py);
101  double cos_phi_x = dx / sqrt(dx * dx + dy * dy);
102  double sin_phi_x = dy / sqrt(dx * dx + dy * dy);
103 
104  double sin_theta_p = pz / sqrt(px * px + py * py + pz * pz);
105  double sin_theta_x = dz / sqrt(dx * dx + dy * dy + dz * dz);
106 
107  double cos_theta_p = sqrt(px * px + py * py) / sqrt(px * px + py * py + pz * pz);
108  double cos_theta_x = sqrt(dx * dx + dy * dy) / sqrt(dx * dx + dy * dy + dz * dz);
109 
110  float feq = sin_phi_p * cos_phi_x - cos_phi_p * sin_phi_x;
111  float seq = sin_theta_p * cos_theta_x - cos_theta_p * sin_theta_x;
112 
113  vl(1) = feq;
114  vl(2) = seq;
115 
116  return std::pair<AlgebraicVector, AlgebraicVector>(vl, point);
117 }
118 
119 std::pair<AlgebraicMatrix, AlgebraicVector> SmartPointingConstraint::makeDerivative(
120  const AlgebraicVector& exPoint) const {
121  AlgebraicMatrix dr(2, 7, 0);
122  AlgebraicVector point = exPoint;
123  double dx = point(1) - refPoint.x();
124  double dy = point(2) - refPoint.y();
125  double dz = point(3) - refPoint.z();
126  double px = point(4);
127  double py = point(5);
128  double pz = point(6);
129 
130  //angular functuions:
131 
132  //half angle solution
133  //d/dx_i
134  dr(1, 1) = (dy * (dx * px + dy * py)) / (pow(pow(dx, 2) + pow(dy, 2), 1.5) * sqrt(pow(px, 2) + pow(py, 2)));
135 
136  dr(1, 2) = -((dx * (dx * px + dy * py)) / (pow(pow(dx, 2) + pow(dy, 2), 1.5) * sqrt(pow(px, 2) + pow(py, 2))));
137 
138  dr(1, 3) = 0;
139 
140  //d/dp_i
141  //debug: x->p index xhange in denominator
142  dr(1, 4) = -((py * (dx * px + dy * py)) / (sqrt(pow(dx, 2) + pow(dy, 2)) * pow(pow(px, 2) + pow(py, 2), 1.5)));
143 
144  dr(1, 5) = (px * (dx * px + dy * py)) / (sqrt(pow(dx, 2) + pow(dy, 2)) * pow(pow(px, 2) + pow(py, 2), 1.5));
145 
146  dr(1, 6) = 0;
147  dr(1, 7) = 0;
148 
149  //2nd equation
150  //d/dx_i
151 
152  dr(2, 1) = (dx * dz * (sqrt(pow(dx, 2) + pow(dy, 2)) * sqrt(pow(px, 2) + pow(py, 2)) + dz * pz)) /
153  (sqrt(pow(dx, 2) + pow(dy, 2)) * pow(pow(dx, 2) + pow(dy, 2) + pow(dz, 2), 1.5) *
154  sqrt(pow(px, 2) + pow(py, 2) + pow(pz, 2)));
155 
156  dr(2, 2) = (dy * dz * (sqrt(pow(dx, 2) + pow(dy, 2)) * sqrt(pow(px, 2) + pow(py, 2)) + dz * pz)) /
157  (sqrt(pow(dx, 2) + pow(dy, 2)) * pow(pow(dx, 2) + pow(dy, 2) + pow(dz, 2), 1.5) *
158  sqrt(pow(px, 2) + pow(py, 2) + pow(pz, 2)));
159 
160  dr(2, 3) = (-((pow(dx, 2) + pow(dy, 2)) * sqrt(pow(px, 2) + pow(py, 2))) - sqrt(pow(dx, 2) + pow(dy, 2)) * dz * pz) /
161  (pow(pow(dx, 2) + pow(dy, 2) + pow(dz, 2), 1.5) * sqrt(pow(px, 2) + pow(py, 2) + pow(pz, 2)));
162 
163  //d/dp_i
164  //debug: x->p index xhange in denominator
165 
166  dr(2, 4) = -((px * pz * (sqrt(pow(dx, 2) + pow(dy, 2)) * sqrt(pow(px, 2) + pow(py, 2)) + dz * pz)) /
167  (sqrt(pow(dx, 2) + pow(dy, 2) + pow(dz, 2)) * sqrt(pow(px, 2) + pow(py, 2)) *
168  pow(pow(px, 2) + pow(py, 2) + pow(pz, 2), 1.5)));
169 
170  dr(2, 5) = -((py * pz * (sqrt(pow(dx, 2) + pow(dy, 2)) * sqrt(pow(px, 2) + pow(py, 2)) + dz * pz)) /
171  (sqrt(pow(dx, 2) + pow(dy, 2) + pow(dz, 2)) * sqrt(pow(px, 2) + pow(py, 2)) *
172  pow(pow(px, 2) + pow(py, 2) + pow(pz, 2), 1.5)));
173 
174  dr(2, 6) = (sqrt(pow(dx, 2) + pow(dy, 2)) * (pow(px, 2) + pow(py, 2)) + dz * sqrt(pow(px, 2) + pow(py, 2)) * pz) /
175  (sqrt(pow(dx, 2) + pow(dy, 2) + pow(dz, 2)) * pow(pow(px, 2) + pow(py, 2) + pow(pz, 2), 1.5));
176 
177  dr(2, 7) = 0;
178 
179  // cout<<"derivative matrix "<<dr<<endl;
180  return std::pair<AlgebraicMatrix, AlgebraicVector>(dr, point);
181 }
SmartPointingConstraint::numberOfEquations
int numberOfEquations() const override
Definition: SmartPointingConstraint.cc:84
cmsswSequenceInfo.seq
seq
Definition: cmsswSequenceInfo.py:539
VertexException
Common base class.
Definition: VertexException.h:12
PV3DBase::x
T x() const
Definition: PV3DBase.h:59
multPhiCorr_741_25nsDY_cfi.py
py
Definition: multPhiCorr_741_25nsDY_cfi.py:12
SmartPointingConstraint::derivative
std::pair< AlgebraicMatrix, AlgebraicVector > derivative(const AlgebraicVector &exPoint) const override
Definition: SmartPointingConstraint.cc:26
PV3DBase::z
T z() const
Definition: PV3DBase.h:61
AlgebraicVector
CLHEP::HepVector AlgebraicVector
Definition: AlgebraicObjects.h:13
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
SmartPointingConstraint::deviations
AlgebraicVector deviations(int nStates) const override
Definition: SmartPointingConstraint.cc:82
SmartPointingConstraint::makeDerivative
std::pair< AlgebraicMatrix, AlgebraicVector > makeDerivative(const AlgebraicVector &exPoint) const
Definition: SmartPointingConstraint.cc:119
PV3DBase::y
T y() const
Definition: PV3DBase.h:60
SmartPointingConstraint::makeValue
std::pair< AlgebraicVector, AlgebraicVector > makeValue(const AlgebraicVector &exPoint) const
Definition: SmartPointingConstraint.cc:86
PVValHelper::dy
Definition: PVValidationHelpers.h:49
multPhiCorr_741_25nsDY_cfi.px
px
Definition: multPhiCorr_741_25nsDY_cfi.py:10
AlgebraicMatrix
CLHEP::HepMatrix AlgebraicMatrix
Definition: AlgebraicObjects.h:14
PVValHelper::dz
Definition: PVValidationHelpers.h:50
flavorHistoryFilter_cfi.dr
dr
Definition: flavorHistoryFilter_cfi.py:37
VertexException.h
SmartPointingConstraint::refPoint
GlobalPoint refPoint
Definition: SmartPointingConstraint.h:55
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
point
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
Definition: invegas.h:5
SmartPointingConstraint.h
PVValHelper::dx
Definition: PVValidationHelpers.h:48
SmartPointingConstraint::value
std::pair< AlgebraicVector, AlgebraicVector > value(const AlgebraicVector &exPoint) const override
Definition: SmartPointingConstraint.cc:4
Generated for CMSSW Reference Manual by doxygen 1.8.16