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