| 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 methods of th 6 // This is the implementation of methods of the HepRotationY class which
7 // were introduced when ZOOM PhysicsVectors wa 7 // were introduced when ZOOM PhysicsVectors was merged in.
8 // 8 //
9 9
>> 10 #ifdef GNUPRAGMA
>> 11 #pragma implementation
>> 12 #endif
>> 13
10 #include "CLHEP/Vector/RotationY.h" 14 #include "CLHEP/Vector/RotationY.h"
11 #include "CLHEP/Vector/AxisAngle.h" 15 #include "CLHEP/Vector/AxisAngle.h"
12 #include "CLHEP/Vector/EulerAngles.h" 16 #include "CLHEP/Vector/EulerAngles.h"
13 #include "CLHEP/Vector/LorentzRotation.h" 17 #include "CLHEP/Vector/LorentzRotation.h"
14 #include "CLHEP/Units/PhysicalConstants.h" 18 #include "CLHEP/Units/PhysicalConstants.h"
15 19
16 #include <cmath> 20 #include <cmath>
17 #include <stdlib.h> 21 #include <stdlib.h>
18 #include <iostream> 22 #include <iostream>
19 23
20 namespace CLHEP { 24 namespace CLHEP {
21 25
22 static inline double safe_acos (double x) { 26 static inline double safe_acos (double x) {
23 if (std::abs(x) <= 1.0) return std::acos(x); 27 if (std::abs(x) <= 1.0) return std::acos(x);
24 return ( (x>0) ? 0 : CLHEP::pi ); 28 return ( (x>0) ? 0 : CLHEP::pi );
25 } 29 }
26 30
27 HepRotationY::HepRotationY(double ddelta) : 31 HepRotationY::HepRotationY(double ddelta) :
28 its_d(proper(ddelta)), its_s(std::sin(ddel 32 its_d(proper(ddelta)), its_s(std::sin(ddelta)), its_c(std::cos(ddelta))
29 {} 33 {}
30 34
31 HepRotationY & HepRotationY::set ( double ddel 35 HepRotationY & HepRotationY::set ( double ddelta ) {
32 its_d = proper(ddelta); 36 its_d = proper(ddelta);
33 its_s = std::sin(its_d); 37 its_s = std::sin(its_d);
34 its_c = std::cos(its_d); 38 its_c = std::cos(its_d);
35 return *this; 39 return *this;
36 } 40 }
37 41
38 double HepRotationY::phi() const { 42 double HepRotationY::phi() const {
39 if ( its_d == 0 ) { 43 if ( its_d == 0 ) {
40 return 0; 44 return 0;
41 } else if ( (its_d < 0) || (its_d == CLHEP:: 45 } else if ( (its_d < 0) || (its_d == CLHEP::pi) ) {
42 return +CLHEP::halfpi; 46 return +CLHEP::halfpi;
43 } else { 47 } else {
44 return -CLHEP::halfpi; 48 return -CLHEP::halfpi;
45 } 49 }
46 } // HepRotationY::phi() 50 } // HepRotationY::phi()
47 51
48 double HepRotationY::theta() const { 52 double HepRotationY::theta() const {
49 return std::fabs( its_d ); 53 return std::fabs( its_d );
50 } // HepRotationY::theta() 54 } // HepRotationY::theta()
51 55
52 double HepRotationY::psi() const { 56 double HepRotationY::psi() const {
53 if ( its_d == 0 ) { 57 if ( its_d == 0 ) {
54 return 0; 58 return 0;
55 } else if ( (its_d < 0) || (its_d == CLHEP:: 59 } else if ( (its_d < 0) || (its_d == CLHEP::pi) ) {
56 return -CLHEP::halfpi; 60 return -CLHEP::halfpi;
57 } else { 61 } else {
58 return +CLHEP::halfpi; 62 return +CLHEP::halfpi;
59 } 63 }
60 } // HepRotationY::psi() 64 } // HepRotationY::psi()
61 65
62 HepEulerAngles HepRotationY::eulerAngles() con 66 HepEulerAngles HepRotationY::eulerAngles() const {
63 return HepEulerAngles( phi(), theta(), ps 67 return HepEulerAngles( phi(), theta(), psi() );
64 } // HepRotationY::eulerAngles() 68 } // HepRotationY::eulerAngles()
65 69
66 70
67 // From the defining code in the implementatio 71 // From the defining code in the implementation of CLHEP (in Rotation.cc)
68 // it is clear that thetaX, phiX form the pola 72 // it is clear that thetaX, phiX form the polar angles in the original
69 // coordinate system of the new X axis (and si 73 // coordinate system of the new X axis (and similarly for phiY and phiZ).
70 // 74 //
71 // This code is taken directly from the origin 75 // This code is taken directly from the original CLHEP. However, there are as
72 // shown opportunities for significant speed i 76 // shown opportunities for significant speed improvement.
73 77
74 double HepRotationY::phiX() const { 78 double HepRotationY::phiX() const {
75 return (yx() == 0.0 && xx() == 0.0) ? 0.0 : 79 return (yx() == 0.0 && xx() == 0.0) ? 0.0 : std::atan2(yx(),xx());
76 // or ---- return 0; 80 // or ---- return 0;
77 } 81 }
78 82
79 double HepRotationY::phiY() const { 83 double HepRotationY::phiY() const {
80 return (yy() == 0.0 && xy() == 0.0) ? 0.0 : 84 return (yy() == 0.0 && xy() == 0.0) ? 0.0 : std::atan2(yy(),xy());
81 // or ---- return CLHEP::halfpi; 85 // or ---- return CLHEP::halfpi;
82 } 86 }
83 87
84 double HepRotationY::phiZ() const { 88 double HepRotationY::phiZ() const {
85 return (yz() == 0.0 && xz() == 0.0) ? 0.0 : 89 return (yz() == 0.0 && xz() == 0.0) ? 0.0 : std::atan2(yz(),xz());
86 // or ---- return 0; 90 // or ---- return 0;
87 } 91 }
88 92
89 double HepRotationY::thetaX() const { 93 double HepRotationY::thetaX() const {
90 return safe_acos(zx()); 94 return safe_acos(zx());
91 } 95 }
92 96
93 double HepRotationY::thetaY() const { 97 double HepRotationY::thetaY() const {
94 return safe_acos(zy()); 98 return safe_acos(zy());
95 // or ---- return CLHEP::halfpi; 99 // or ---- return CLHEP::halfpi;
96 } 100 }
97 101
98 double HepRotationY::thetaZ() const { 102 double HepRotationY::thetaZ() const {
99 return safe_acos(zz()); 103 return safe_acos(zz());
100 // or ---- return d; 104 // or ---- return d;
101 } 105 }
102 106
103 void HepRotationY::setDelta ( double ddelta ) 107 void HepRotationY::setDelta ( double ddelta ) {
104 set(ddelta); 108 set(ddelta);
105 } 109 }
106 110
107 void HepRotationY::decompose 111 void HepRotationY::decompose
108 (HepAxisAngle & rotation, Hep3Vector & boost 112 (HepAxisAngle & rotation, Hep3Vector & boost) const {
109 boost.set(0,0,0); 113 boost.set(0,0,0);
110 rotation = axisAngle(); 114 rotation = axisAngle();
111 } 115 }
112 116
113 void HepRotationY::decompose 117 void HepRotationY::decompose
114 (Hep3Vector & boost, HepAxisAngle & rotation 118 (Hep3Vector & boost, HepAxisAngle & rotation) const {
115 boost.set(0,0,0); 119 boost.set(0,0,0);
116 rotation = axisAngle(); 120 rotation = axisAngle();
117 } 121 }
118 122
119 void HepRotationY::decompose 123 void HepRotationY::decompose
120 (HepRotation & rotation, HepBoost & bo 124 (HepRotation & rotation, HepBoost & boost) const {
121 boost.set(0,0,0); 125 boost.set(0,0,0);
122 rotation = HepRotation(*this); 126 rotation = HepRotation(*this);
123 } 127 }
124 128
125 void HepRotationY::decompose 129 void HepRotationY::decompose
126 (HepBoost & boost, HepRotation & rotat 130 (HepBoost & boost, HepRotation & rotation) const {
127 boost.set(0,0,0); 131 boost.set(0,0,0);
128 rotation = HepRotation(*this); 132 rotation = HepRotation(*this);
129 } 133 }
130 134
131 double HepRotationY::distance2( const HepRotat 135 double HepRotationY::distance2( const HepRotationY & r ) const {
132 double answer = 2.0 * ( 1.0 - ( its_s * r.it 136 double answer = 2.0 * ( 1.0 - ( its_s * r.its_s + its_c * r.its_c ) ) ;
133 return (answer >= 0) ? answer : 0; 137 return (answer >= 0) ? answer : 0;
134 } 138 }
135 139
136 double HepRotationY::distance2( const HepRotat 140 double HepRotationY::distance2( const HepRotation & r ) const {
137 double sum = xx() * r.xx() + 141 double sum = xx() * r.xx() + xz() * r.xz()
138 + r.yy() 142 + r.yy()
139 + zx() * r.zx() 143 + zx() * r.zx() + zz() * r.zz();
140 double answer = 3.0 - sum; 144 double answer = 3.0 - sum;
141 return (answer >= 0 ) ? answer : 0; 145 return (answer >= 0 ) ? answer : 0;
142 } 146 }
143 147
144 double HepRotationY::distance2( const HepLoren 148 double HepRotationY::distance2( const HepLorentzRotation & lt ) const {
145 HepAxisAngle a; 149 HepAxisAngle a;
146 Hep3Vector b; 150 Hep3Vector b;
147 lt.decompose(b, a); 151 lt.decompose(b, a);
148 double bet = b.beta(); 152 double bet = b.beta();
149 double bet2 = bet*bet; 153 double bet2 = bet*bet;
150 HepRotation r(a); 154 HepRotation r(a);
151 return bet2/(1-bet2) + distance2(r); 155 return bet2/(1-bet2) + distance2(r);
152 } 156 }
153 157
154 double HepRotationY::distance2( const HepBoost 158 double HepRotationY::distance2( const HepBoost & lt ) const {
155 return distance2( HepLorentzRotation(lt)); 159 return distance2( HepLorentzRotation(lt));
156 } 160 }
157 161
158 double HepRotationY::howNear( const HepRotatio 162 double HepRotationY::howNear( const HepRotationY & r ) const {
159 return std::sqrt(distance2(r)); 163 return std::sqrt(distance2(r));
160 } 164 }
161 double HepRotationY::howNear( const HepRotatio 165 double HepRotationY::howNear( const HepRotation & r ) const {
162 return std::sqrt(distance2(r)); 166 return std::sqrt(distance2(r));
163 } 167 }
164 double HepRotationY::howNear( const HepBoost & 168 double HepRotationY::howNear( const HepBoost & lt ) const {
165 return std::sqrt(distance2(lt)); 169 return std::sqrt(distance2(lt));
166 } 170 }
167 double HepRotationY::howNear( const HepLorentz 171 double HepRotationY::howNear( const HepLorentzRotation & lt ) const {
168 return std::sqrt(distance2(lt)); 172 return std::sqrt(distance2(lt));
169 } 173 }
170 bool HepRotationY::isNear(const HepRotationY & 174 bool HepRotationY::isNear(const HepRotationY & r,double epsilon)const{
171 return (distance2(r) <= epsilon*epsilon); 175 return (distance2(r) <= epsilon*epsilon);
172 } 176 }
173 bool HepRotationY::isNear(const HepRotation & 177 bool HepRotationY::isNear(const HepRotation & r,double epsilon)const {
174 return (distance2(r) <= epsilon*epsilon); 178 return (distance2(r) <= epsilon*epsilon);
175 } 179 }
176 bool HepRotationY::isNear( const HepBoost & lt 180 bool HepRotationY::isNear( const HepBoost & lt,double epsilon) const {
177 return (distance2(lt) <= epsilon*epsilon); 181 return (distance2(lt) <= epsilon*epsilon);
178 } 182 }
179 bool HepRotationY::isNear( const HepLorentzRot 183 bool HepRotationY::isNear( const HepLorentzRotation & lt,
180 double ep 184 double epsilon) const {
181 return (distance2(lt) <= epsilon*epsilon); 185 return (distance2(lt) <= epsilon*epsilon);
182 } 186 }
183 187
184 double HepRotationY::norm2() const { 188 double HepRotationY::norm2() const {
185 return 2.0 - 2.0 * its_c; 189 return 2.0 - 2.0 * its_c;
186 } 190 }
187 191
188 std::ostream & HepRotationY::print( std::ostre 192 std::ostream & HepRotationY::print( std::ostream & os ) const {
189 os << "\nRotation about Y (" << its_d << 193 os << "\nRotation about Y (" << its_d <<
190 ") [cos d = " << its_c << " si 194 ") [cos d = " << its_c << " sin d = " << its_s << "]\n";
191 return os; 195 return os;
192 } 196 }
193 197
194 } // namespace CLHEP 198 } // namespace CLHEP
195 199