LCFIVertex  0.7.2
vector3.cpp
1 #include "../inc/vector3.h"
2 #include "../inc/projection.h"
3 #include <boost/version.hpp>
4 
5 namespace vertex_lcfi
6 {
7 namespace util
8 {
9 
10  //Default Constructor
11  Vector3::Vector3()
12  {
13  Base_Vector::clear();
14  }
15 
16  void Vector3::operator=(const Base_Vector &r)
17  {
18  Base_Vector::operator=(r);
19  }
20 
21  bool Vector3::operator!=(const Vector3 &r)
22  {
23  return (r.x()!=x() || r.y()!=y() || r.z()!=z());
24  }
25 
26  // Construction with components
27  Vector3::Vector3(const double &c1,const double &c2,const double &c3)
28  {
29 #if BOOST_VERSION>=103300
30  insert_element(0,c1);
31  insert_element(1,c2);
32  insert_element(2,c3);
33 #else
34  insert(0,c1);
35  insert(1,c2);
36  insert(2,c3);
37 #endif
38  }
39  // Accessor functions
40  double Vector3::x() const
41  {
42  return operator()(0);
43  }
44  double Vector3::y() const
45  {
46  return operator()(1);
47  }
48  double Vector3::z() const
49  {
50  return operator()(2);
51  }
52  double & Vector3::x()
53  {
54  return operator()(0);
55  }
56  double & Vector3::y()
57  {
58  return operator()(1);
59  }
60  double & Vector3::z()
61  {
62  return operator()(2);
63  }
64  double Vector3::px() const
65  {
66  return operator()(0);
67  }
68  double Vector3::py() const
69  {
70  return operator()(1);
71  }
72  double Vector3::pz() const
73  {
74  return operator()(2);
75  }
76  double & Vector3::px()
77  {
78  return operator()(0);
79  }
80  double & Vector3::py()
81  {
82  return operator()(1);
83  }
84  double & Vector3::pz()
85  {
86  return operator()(2);
87  }
88  double Vector3::mag2() const
89  {
90  return boost::numeric::ublas::inner_prod(*this,*this);
91  }
92  double Vector3::mag() const
93  {
94  return std::sqrt(boost::numeric::ublas::inner_prod(*this,*this));
95  }
96  double Vector3::mag(Projection Proj) const
97  {
98  switch (Proj)
99  {
100  case ThreeD:
101  return std::sqrt(boost::numeric::ublas::inner_prod(*this,*this));
102  case RPhi:
103  return std::sqrt(pow(this->x(),2) + pow(this->y(),2));
104  case Z:
105  return fabs(this->z());
106  }
107  //TODO Throw Something
108  return 0;
109  }
110  double Vector3::mag2(Projection Proj) const
111  {
112  switch (Proj)
113  {
114  case ThreeD:
115  return boost::numeric::ublas::inner_prod(*this,*this);
116  case RPhi:
117  return pow(this->x(),2) + pow(this->y(),2);
118  case Z:
119  return pow(this->z(),2);
120  }
121  //TODO Throw Something
122  return 0;
123  }
124 
125  double Vector3::p2() const
126  {
127  return mag2();
128  }
129  double Vector3::p() const
130  {
131  return mag();
132  }
133 
134  double Vector3::r() const
135  {
136  return mag();
137  }
138  double Vector3::theta() const
139  {
140  return operator()(0) == 0.0 && operator()(1) == 0.0 && operator()(2) == 0.0 ? 0.0 : std::atan2(perp(),operator()(2));
141  }
142  double Vector3::phi() const
143  {
144  return operator()(0) == 0.0 && operator()(1) == 0.0 ? 0.0 : std::atan2(operator()(1),operator()(0));
145  }
146  double Vector3::cosTheta() const
147  {
148  const double ptot = mag();
149  return ptot == 0.0 ? 1.0 : operator()(2)/ptot;
150  }
151 
152  double Vector3::perp2() const
153  {
154  return operator()(0)*operator()(0)+operator()(1)*operator()(1);
155  }
156  double Vector3::perp() const
157  {
158  return std::sqrt(perp2());
159  }
160 
161  double Vector3::pt2() const
162  {
163  return perp2();
164  }
165  double Vector3::pt() const
166  {
167  return perp();
168  }
169 
170  double Vector3::eta() const
171  {
172  const double m = mag();
173  if (m == 0.0) return 0.0;
174  const double tz = z();
175  if (m == tz) return 1.0E72;
176  if (m == -tz) return -1.0E72;
177  return 0.5*std::log( (m+tz)/(m-tz) );
178  }
179 
180  // Turn this vector into a unit vector with the same direction
181  void Vector3::makeUnit()
182  {
183  const double m = mag();
184  if (m != 0.0)
185  {
186  x() /= m;
187  y() /= m;
188  z() /= m;
189  }
190  }
191 
192  Vector3 Vector3::unit() const
193  {
194  const double m = mag();
195  if (m != 0.0)
196  return Vector3(x() / m,y() / m,z() / m);
197  else
198  return Vector3(x(),y(),z());
199  }
200 
201  // Dot product with another Vector3
202  double Vector3::dot(const Vector3 &v) const
203  {
204  return boost::numeric::ublas::inner_prod(*this,v);
205  }
206 
207  // Cross product with another Vector3. Result is (this x v)
208  Vector3 Vector3::cross(const Vector3 &v) const
209  {
210  return Vector3(y()*v.z()-z()*v.y(),z()*v.x()-x()*v.z(),x()*v.y()-y()*v.x());
211  }
212 
213  Vector3 Vector3::subtract(const Vector3 &v) const
214  {
215  return Vector3(x()-v.x(), y()-v.y(), z()-v.z());
216  }
217 
218  Vector3 Vector3::add(const Vector3 &v) const
219  {
220  return Vector3(x()+v.x(), y()+v.y(), z()+v.z());
221  }
222 
223  Vector3 Vector3::mult(double k) const
224  {
225  return Vector3(x()*k, y()*k, z()*k);
226  }
227 
228  double Vector3::distanceTo(const Vector3 &v) const
229  {
230  return (this->subtract(v)).mag();
231  }
232 
233  double Vector3::distanceTo2(const Vector3 &v) const
234  {
235  return (this->subtract(v)).mag2();
236  }
237 
238  double Vector3::distanceToLine(const Vector3 &Point, const Vector3 &Direction) const
239  {
240  return std::sqrt(distanceToLine2(Point, Direction));
241  }
242 
243  double Vector3::distanceToLine2(const Vector3 &Point, const Vector3 &Direction) const
244  {
245  if (Direction.mag()==0)
246  return -1;
247  //TODO Exception above??
248  Vector3 P1 = Point;
249  Vector3 P2 = Point+Direction;
250  Vector3 P3 = Vector3(x(),y(),z());
251 
252  //Parametric distance along line to perp to point
253  double u = (((P3.x()-P1.x())*(P2.x()-P1.x())) + ((P3.y()-P1.y())*(P2.y()-P1.y())) + ((P3.z()-P1.z())*(P2.z()-P1.z()))) / (P2.subtract(P1)).mag2();
254  Vector3 PerpPoint = P1+((P2-P1)*u);
255 
256  return PerpPoint.distanceTo2(P3);
257  }
258 
259  std::vector<double> Vector3::stlVector() const
260  {
261  std::vector<double> out;
262  out.push_back(x());
263  out.push_back(y());
264  out.push_back(z());
265  return out;
266  }
267 
268 
269 
270 }
271 }