CLHEP VERSION Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // -*- C++ -*- 00002 // --------------------------------------------------------------------------- 00003 // 00004 // This file is a part of the CLHEP - a Class Library for High Energy Physics. 00005 // 00006 // This is the implementation of the Hep2Vector class. 00007 // 00008 //------------------------------------------------------------- 00009 00010 #include "CLHEP/Vector/defs.h" 00011 #include "CLHEP/Vector/TwoVector.h" 00012 #include "CLHEP/Vector/ZMxpv.h" 00013 #include "CLHEP/Vector/ThreeVector.h" 00014 00015 #include <cmath> 00016 #include <iostream> 00017 00018 namespace CLHEP { 00019 00020 double Hep2Vector::tolerance = Hep2Vector::ZMpvToleranceTicks * 2.22045e-16; 00021 00022 double Hep2Vector::setTolerance (double tol) { 00023 // Set the tolerance for Hep2Vectors to be considered near one another 00024 double oldTolerance (tolerance); 00025 tolerance = tol; 00026 return oldTolerance; 00027 } 00028 00029 double Hep2Vector::operator () (int i) const { 00030 if (i == 0) { 00031 return x(); 00032 }else if (i == 1) { 00033 return y(); 00034 }else{ 00035 ZMthrowA(ZMxpvIndexRange( 00036 "Hep2Vector::operator(): bad index")); 00037 return 0.0; 00038 } 00039 } 00040 00041 double & Hep2Vector::operator () (int i) { 00042 static double dummy; 00043 switch(i) { 00044 case X: 00045 return dx; 00046 case Y: 00047 return dy; 00048 default: 00049 ZMthrowA (ZMxpvIndexRange( 00050 "Hep2Vector::operator() : bad index")); 00051 return dummy; 00052 } 00053 } 00054 00055 void Hep2Vector::rotate(double angler) { 00056 double s1 = std::sin(angler); 00057 double c = std::cos(angler); 00058 double xx = dx; 00059 dx = c*xx - s1*dy; 00060 dy = s1*xx + c*dy; 00061 } 00062 00063 Hep2Vector operator/ (const Hep2Vector & p, double a) { 00064 if (a==0) { 00065 ZMthrowA(ZMxpvInfiniteVector( "Division of Hep2Vector by zero")); 00066 } 00067 return Hep2Vector(p.x()/a, p.y()/a); 00068 } 00069 00070 std::ostream & operator << (std::ostream & os, const Hep2Vector & q) { 00071 os << "(" << q.x() << ", " << q.y() << ")"; 00072 return os; 00073 } 00074 00075 void ZMinput2doubles ( std::istream & is, const char * type, 00076 double & x, double & y ); 00077 00078 std::istream & operator>>(std::istream & is, Hep2Vector & p) { 00079 double x, y; 00080 ZMinput2doubles ( is, "Hep2Vector", x, y ); 00081 p.set(x, y); 00082 return is; 00083 } // operator>>() 00084 00085 Hep2Vector::operator Hep3Vector () const { 00086 return Hep3Vector ( dx, dy, 0.0 ); 00087 } 00088 00089 int Hep2Vector::compare (const Hep2Vector & v) const { 00090 if ( dy > v.dy ) { 00091 return 1; 00092 } else if ( dy < v.dy ) { 00093 return -1; 00094 } else if ( dx > v.dx ) { 00095 return 1; 00096 } else if ( dx < v.dx ) { 00097 return -1; 00098 } else { 00099 return 0; 00100 } 00101 } /* Compare */ 00102 00103 00104 bool Hep2Vector::operator > (const Hep2Vector & v) const { 00105 return (compare(v) > 0); 00106 } 00107 bool Hep2Vector::operator < (const Hep2Vector & v) const { 00108 return (compare(v) < 0); 00109 } 00110 bool Hep2Vector::operator>= (const Hep2Vector & v) const { 00111 return (compare(v) >= 0); 00112 } 00113 bool Hep2Vector::operator<= (const Hep2Vector & v) const { 00114 return (compare(v) <= 0); 00115 } 00116 00117 bool Hep2Vector::isNear(const Hep2Vector & p, double epsilon) const { 00118 double limit = dot(p)*epsilon*epsilon; 00119 return ( (*this - p).mag2() <= limit ); 00120 } /* isNear() */ 00121 00122 double Hep2Vector::howNear(const Hep2Vector & p ) const { 00123 double d = (*this - p).mag2(); 00124 double pdp = dot(p); 00125 if ( (pdp > 0) && (d < pdp) ) { 00126 return std::sqrt (d/pdp); 00127 } else if ( (pdp == 0) && (d == 0) ) { 00128 return 0; 00129 } else { 00130 return 1; 00131 } 00132 } /* howNear */ 00133 00134 double Hep2Vector::howParallel (const Hep2Vector & v) const { 00135 // | V1 x V2 | / | V1 dot V2 | 00136 // Of course, the "cross product" is fictitious but the math is valid 00137 double v1v2 = std::fabs(dot(v)); 00138 if ( v1v2 == 0 ) { 00139 // Zero is parallel to no other vector except for zero. 00140 return ( (mag2() == 0) && (v.mag2() == 0) ) ? 0 : 1; 00141 } 00142 double abscross = std::fabs ( dx * v.y() - dy - v.x() ); 00143 if ( abscross >= v1v2 ) { 00144 return 1; 00145 } else { 00146 return abscross/v1v2; 00147 } 00148 } /* howParallel() */ 00149 00150 bool Hep2Vector::isParallel (const Hep2Vector & v, 00151 double epsilon) const { 00152 // | V1 x V2 | <= epsilon * | V1 dot V2 | 00153 // Of course, the "cross product" is fictitious but the math is valid 00154 double v1v2 = std::fabs(dot(v)); 00155 if ( v1v2 == 0 ) { 00156 // Zero is parallel to no other vector except for zero. 00157 return ( (mag2() == 0) && (v.mag2() == 0) ); 00158 } 00159 double abscross = std::fabs ( dx * v.y() - dy - v.x() ); 00160 return ( abscross <= epsilon * v1v2 ); 00161 } /* isParallel() */ 00162 00163 double Hep2Vector::howOrthogonal (const Hep2Vector & v) const { 00164 // | V1 dot V2 | / | V1 x V2 | 00165 // Of course, the "cross product" is fictitious but the math is valid 00166 double v1v2 = std::fabs(dot(v)); 00167 if ( v1v2 == 0 ) { 00168 return 0; // Even if one or both are 0, they are considered orthogonal 00169 } 00170 double abscross = std::fabs ( dx * v.y() - dy - v.x() ); 00171 if ( v1v2 >= abscross ) { 00172 return 1; 00173 } else { 00174 return v1v2/abscross; 00175 } 00176 } /* howOrthogonal() */ 00177 00178 bool Hep2Vector::isOrthogonal (const Hep2Vector & v, 00179 double epsilon) const { 00180 // | V1 dot V2 | <= epsilon * | V1 x V2 | 00181 // Of course, the "cross product" is fictitious but the math is valid 00182 double v1v2 = std::fabs(dot(v)); 00183 double abscross = std::fabs ( dx * v.y() - dy - v.x() ); 00184 return ( v1v2 <= epsilon * abscross ); 00185 } /* isOrthogonal() */ 00186 00187 } // namespace CLHEP