CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ColinearityKinematicConstraint.cc
Go to the documentation of this file.
3 
5 {
6  dimension = dim;
7  if (dimension == Phi) size = 1;
8  else size = 2;
9 }
10 
12  const GlobalPoint& point) const
13 {
14  if(states.size()<2) throw VertexException("ColinearityKinematicConstraint::<2 states passed");
15  AlgebraicVector res(size,0);
16 
17  double a_1 = -states[0].particleCharge()*states[0].magneticField()->inInverseGeV(states[0].globalPosition()).z();
18  double a_2 = -states[1].particleCharge()*states[1].magneticField()->inInverseGeV(states[1].globalPosition()).z();
19 
20  AlgebraicVector7 p1 = states[0].kinematicParameters().vector();
21  AlgebraicVector7 p2 = states[1].kinematicParameters().vector();
22 
23  double p1vx = p1(3) - a_1*(point.y() - p1(1));
24  double p1vy = p1(4) + a_1*(point.x() - p1(0));
25  double pt1 = sqrt(p1(3)*p1(3)+p1(4)*p1(4));
26 
27  double p2vx = p2(3) - a_2*(point.y() - p2(1));
28  double p2vy = p2(4) + a_2*(point.x() - p2(0));
29  double pt2 = sqrt(p2(3)*p2(3)+p2(4)*p2(4));
30 
31  // H_phi:
32  res(1) = atan2(p1vy,p1vx) - atan2(p2vy,p2vx);
33  if ( res(1) > M_PI ) res(1) -= 2.0*M_PI;
34  if ( res(1) <= -M_PI ) res(1) += 2.0*M_PI;
35  // H_theta:
36  if (dimension == PhiTheta) {
37  res(2) = atan2(pt1,p1(5)) - atan2(pt2,p2(5));
38  if ( res(2) > M_PI ) res(2) -= 2.0*M_PI;
39  if ( res(2) <= -M_PI ) res(2) += 2.0*M_PI;
40  }
41 
42 // cout << res(1) << " "<<res(2) << "\n ";// res(2)
43 
44  return res;
45 }
46 
48  const GlobalPoint& point) const
49 {
50  int n_st = states.size();
51  if(n_st<2) throw VertexException("ColinearityKinematicConstraint::<2 states passed");
52  AlgebraicMatrix res(size,n_st*7,0);
53 
54  double a_1 = -states[0].particleCharge()*states[0].magneticField()->inInverseGeV(states[0].globalPosition()).z();
55  double a_2 = -states[1].particleCharge()*states[1].magneticField()->inInverseGeV(states[1].globalPosition()).z();
56 
57  AlgebraicVector7 p1 = states[0].kinematicParameters().vector();
58  AlgebraicVector7 p2 = states[1].kinematicParameters().vector();
59 
60  double p1vx = p1(3) - a_1*(point.y() - p1(1));
61  double p1vy = p1(4) + a_1*(point.x() - p1(0));
62  double k1 = 1.0/(p1vx*p1vx + p1vy*p1vy);
63  double pt1 = sqrt(p1(3)*p1(3)+p1(4)*p1(4));
64  double pTot1 = sqrt(p1(3)*p1(3)+p1(4)*p1(4)+p1(5)*p1(5));
65 
66  double p2vx = p2(3) - a_2*(point.y() - p2(1));
67  double p2vy = p2(4) + a_2*(point.x() - p2(0));
68  double k2 = 1.0/(p2vx*p2vx + p2vy*p2vy);
69  double pt2 = sqrt(p2(3)*p2(3)+p2(4)*p2(4));
70  double pTot2 = sqrt(p2(3)*p2(3)+p2(4)*p2(4)+p2(5)*p2(5));
71 
72  // H_phi:
73 
74  //x1 and x2 derivatives: 1st and 8th elements
75  res(1,1) = -k1*a_1*p1vx;
76  res(1,8) = k2*a_2*p2vx;
77 
78  //y1 and y2 derivatives: 2nd and 9th elements:
79  res(1,2) = -k1*a_1*p1vy;
80  res(1,9) = k2*a_2*p2vy;
81 
82  //z1 and z2 components: 3d and 10th elmnets stay 0:
83  res(1,3) = 0.; res(1,10) = 0.;
84 
85  //px1 and px2 components: 4th and 11th elements:
86  res(1,4) = -k1*p1vy;
87  res(1,11) = k2*p2vy;
88 
89  //py1 and py2 components: 5th and 12 elements:
90  res(1,5) = k1*p1vx;
91  res(1,12) = -k2*p2vx;
92 
93 
94  //pz1 and pz2 components: 6th and 13 elements:
95  res(1,6) = 0.; res(1,13) = 0.;
96  //mass components: 7th and 14th elements:
97  res(1,7) = 0.; res(1,14) = 0.;
98 
99  if (dimension == PhiTheta) {
100  // H_theta:
101  //x1 and x2 derivatives: 1st and 8th elements
102  res(2,1) = 0.; res(2,8) = 0.;
103 
104  //y1 and y2 derivatives: 2nd and 9th elements:
105  res(2,2) = 0.; res(2,9) = 0.;
106 
107  //z1 and z2 components: 3d and 10th elmnets stay 0:
108  res(2,3) = 0.; res(2,10) = 0.;
109 
110  res(2,4) = p1(5)*p1(3) / (pTot1*pTot1*pt1);
111  res(2,11) = - p2(5)*p2(3) / (pTot2*pTot2*pt2);
112 
113  //py1 and py2 components: 5th and 12 elements:
114  res(2,5) = p1(5)*p1(4) / (pTot1*pTot1*pt1);
115  res(2,12) = - p2(5)*p2(4) / (pTot2*pTot2*pt2);
116 
117  //pz1 and pz2 components: 6th and 13 elements:
118  res(2,6) = - pt1/(pTot1*pTot1);
119  res(2,13) = pt2/(pTot2*pTot2);
120  //mass components: 7th and 14th elements:
121  res(2,7) = 0.; res(2,14) = 0.;
122  }
123 
124  return res;
125 }
126 
128  const GlobalPoint& point) const
129 {
130  AlgebraicMatrix res(size,3,0);
131  if(states.size()<2) throw VertexException("ColinearityKinematicConstraint::<2 states passed");
132 
133  double a_1 = -states[0].particleCharge() * states[0].magneticField()->inInverseGeV(states[0].globalPosition()).z();
134  double a_2 = -states[1].particleCharge() * states[1].magneticField()->inInverseGeV(states[1].globalPosition()).z();
135 
136  AlgebraicVector7 p1 = states[0].kinematicParameters().vector();
137  AlgebraicVector7 p2 = states[1].kinematicParameters().vector();
138 
139  double p1vx = p1(3) - a_1*(point.y() - p1(1));
140  double p1vy = p1(4) + a_1*(point.x() - p1(0));
141  double k1 = 1.0/(p1vx*p1vx + p1vy*p1vy);
142  //double pt1 = sqrt(p1(3)*p1(3)+p1(4)*p1(4));
143 
144  double p2vx = p2(3) - a_2*(point.y() - p2(1));
145  double p2vy = p2(4) + a_2*(point.x() - p2(0));
146  double k2 = 1.0/(p2vx*p2vx + p2vy*p2vy);
147  //double pt2 = sqrt(p2(3)*p2(3)+p2(4)*p2(4));
148 
149  // H_phi:
150 
151  // xv component
152  res(1,1) = k1*a_1*p1vx - k2*a_2*p2vx;
153 
154  //yv component
155  res(1,2) = k1*a_1*p1vy - k2*a_2*p2vy;
156 
157  //zv component
158  res(1,3) = 0.;
159 
160  // H_theta: no correlation with vertex position
161  if (dimension == PhiTheta) {
162  res(2,1) = 0.;
163  res(2,2) = 0.;
164  res(2,3) = 0.;
165  }
166 
167  return res;
168 }
Common base class.
virtual AlgebraicMatrix positionDerivative(const std::vector< KinematicState > &states, const GlobalPoint &point) const
T y() const
Definition: PV3DBase.h:63
ROOT::Math::SVector< double, 7 > AlgebraicVector7
Definition: Matrices.h:15
CLHEP::HepMatrix AlgebraicMatrix
T sqrt(T t)
Definition: SSEVec.h:48
double p2[4]
Definition: TauolaWrapper.h:90
#define M_PI
virtual AlgebraicMatrix parametersDerivative(const std::vector< KinematicState > &states, const GlobalPoint &point) const
CLHEP::HepVector AlgebraicVector
double p1[4]
Definition: TauolaWrapper.h:89
T x() const
Definition: PV3DBase.h:62
*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
virtual AlgebraicVector value(const std::vector< KinematicState > &states, const GlobalPoint &point) const
ColinearityKinematicConstraint(ConstraintDim dim=Phi)