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