CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

TwoVector.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 // ---------------------------------------------------------------------------
3 //
4 // This file is a part of the CLHEP - a Class Library for High Energy Physics.
5 //
6 // This is the implementation of the Hep2Vector class.
7 //
8 //-------------------------------------------------------------
9 
10 #include "CLHEP/Vector/defs.h"
11 #include "CLHEP/Vector/TwoVector.h"
12 #include "CLHEP/Vector/ZMxpv.h"
13 #include "CLHEP/Vector/ThreeVector.h"
14 
15 #include <cmath>
16 #include <iostream>
17 
18 namespace CLHEP {
19 
20 double Hep2Vector::tolerance = Hep2Vector::ZMpvToleranceTicks * 2.22045e-16;
21 
22 double Hep2Vector::setTolerance (double tol) {
23 // Set the tolerance for Hep2Vectors to be considered near one another
24  double oldTolerance (tolerance);
25  tolerance = tol;
26  return oldTolerance;
27 }
28 
29 double Hep2Vector::operator () (int i) const {
30  if (i == 0) {
31  return x();
32  }else if (i == 1) {
33  return y();
34  }else{
35  ZMthrowA(ZMxpvIndexRange(
36  "Hep2Vector::operator(): bad index"));
37  return 0.0;
38  }
39 }
40 
41 double & Hep2Vector::operator () (int i) {
42  static double dummy;
43  switch(i) {
44  case X:
45  return dx;
46  case Y:
47  return dy;
48  default:
49  ZMthrowA (ZMxpvIndexRange(
50  "Hep2Vector::operator() : bad index"));
51  return dummy;
52  }
53 }
54 
55 void Hep2Vector::rotate(double angler) {
56  double s1 = std::sin(angler);
57  double c = std::cos(angler);
58  double xx = dx;
59  dx = c*xx - s1*dy;
60  dy = s1*xx + c*dy;
61 }
62 
63 Hep2Vector operator/ (const Hep2Vector & p, double a) {
64  if (a==0) {
65  ZMthrowA(ZMxpvInfiniteVector( "Division of Hep2Vector by zero"));
66  }
67  return Hep2Vector(p.x()/a, p.y()/a);
68 }
69 
70 std::ostream & operator << (std::ostream & os, const Hep2Vector & q) {
71  os << "(" << q.x() << ", " << q.y() << ")";
72  return os;
73 }
74 
75 void ZMinput2doubles ( std::istream & is, const char * type,
76  double & x, double & y );
77 
78 std::istream & operator>>(std::istream & is, Hep2Vector & p) {
79  double x, y;
80  ZMinput2doubles ( is, "Hep2Vector", x, y );
81  p.set(x, y);
82  return is;
83 } // operator>>()
84 
85 Hep2Vector::operator Hep3Vector () const {
86  return Hep3Vector ( dx, dy, 0.0 );
87 }
88 
89 int Hep2Vector::compare (const Hep2Vector & v) const {
90  if ( dy > v.dy ) {
91  return 1;
92  } else if ( dy < v.dy ) {
93  return -1;
94  } else if ( dx > v.dx ) {
95  return 1;
96  } else if ( dx < v.dx ) {
97  return -1;
98  } else {
99  return 0;
100  }
101 } /* Compare */
102 
103 
104 bool Hep2Vector::operator > (const Hep2Vector & v) const {
105  return (compare(v) > 0);
106 }
107 bool Hep2Vector::operator < (const Hep2Vector & v) const {
108  return (compare(v) < 0);
109 }
110 bool Hep2Vector::operator>= (const Hep2Vector & v) const {
111  return (compare(v) >= 0);
112 }
113 bool Hep2Vector::operator<= (const Hep2Vector & v) const {
114  return (compare(v) <= 0);
115 }
116 
117 bool Hep2Vector::isNear(const Hep2Vector & p, double epsilon) const {
118  double limit = dot(p)*epsilon*epsilon;
119  return ( (*this - p).mag2() <= limit );
120 } /* isNear() */
121 
122 double Hep2Vector::howNear(const Hep2Vector & p ) const {
123  double d = (*this - p).mag2();
124  double pdp = dot(p);
125  if ( (pdp > 0) && (d < pdp) ) {
126  return std::sqrt (d/pdp);
127  } else if ( (pdp == 0) && (d == 0) ) {
128  return 0;
129  } else {
130  return 1;
131  }
132 } /* howNear */
133 
134 double Hep2Vector::howParallel (const Hep2Vector & v) const {
135  // | V1 x V2 | / | V1 dot V2 |
136  // Of course, the "cross product" is fictitious but the math is valid
137  double v1v2 = std::fabs(dot(v));
138  if ( v1v2 == 0 ) {
139  // Zero is parallel to no other vector except for zero.
140  return ( (mag2() == 0) && (v.mag2() == 0) ) ? 0 : 1;
141  }
142  double abscross = std::fabs ( dx * v.y() - dy - v.x() );
143  if ( abscross >= v1v2 ) {
144  return 1;
145  } else {
146  return abscross/v1v2;
147  }
148 } /* howParallel() */
149 
151  double epsilon) const {
152  // | V1 x V2 | <= epsilon * | V1 dot V2 |
153  // Of course, the "cross product" is fictitious but the math is valid
154  double v1v2 = std::fabs(dot(v));
155  if ( v1v2 == 0 ) {
156  // Zero is parallel to no other vector except for zero.
157  return ( (mag2() == 0) && (v.mag2() == 0) );
158  }
159  double abscross = std::fabs ( dx * v.y() - dy - v.x() );
160  return ( abscross <= epsilon * v1v2 );
161 } /* isParallel() */
162 
163 double Hep2Vector::howOrthogonal (const Hep2Vector & v) const {
164  // | V1 dot V2 | / | V1 x V2 |
165  // Of course, the "cross product" is fictitious but the math is valid
166  double v1v2 = std::fabs(dot(v));
167  if ( v1v2 == 0 ) {
168  return 0; // Even if one or both are 0, they are considered orthogonal
169  }
170  double abscross = std::fabs ( dx * v.y() - dy - v.x() );
171  if ( v1v2 >= abscross ) {
172  return 1;
173  } else {
174  return v1v2/abscross;
175  }
176 } /* howOrthogonal() */
177 
179  double epsilon) const {
180  // | V1 dot V2 | <= epsilon * | V1 x V2 |
181  // Of course, the "cross product" is fictitious but the math is valid
182  double v1v2 = std::fabs(dot(v));
183  double abscross = std::fabs ( dx * v.y() - dy - v.x() );
184  return ( v1v2 <= epsilon * abscross );
185 } /* isOrthogonal() */
186 
187 } // namespace CLHEP