14 #pragma implementation
17 #include "CLHEP/Vector/ThreeVector.h"
18 #include "CLHEP/Units/PhysicalConstants.h"
25 void Hep3Vector::setMag(
double ma) {
28 std::cerr <<
"Hep3Vector::setMag() - "
29 <<
"zero vector can't be stretched" << std::endl;
38 Hep3Vector & Hep3Vector::rotateUz(
const Hep3Vector& NewUzVector) {
41 double u1 = NewUzVector.x();
42 double u2 = NewUzVector.y();
43 double u3 = NewUzVector.z();
44 double up = u1*u1 + u2*u2;
48 double px = dx, py = dy, pz = dz;
49 dx = (u1*u3*px - u2*py)/up + u1*pz;
50 dy = (u2*u3*px + u1*py)/up + u2*pz;
53 else if (u3 < 0.) { dx = -dx; dz = -dz; }
58 double Hep3Vector::pseudoRapidity()
const {
60 if ( m1== 0 )
return 0.0;
61 if ( m1==
z() )
return 1.0E72;
62 if ( m1== -
z() )
return -1.0E72;
63 return 0.5*std::log( (m1+
z())/(m1-
z()) );
66 std::ostream &
operator<< (std::ostream & os,
const Hep3Vector & v) {
67 return os <<
"(" << v.x() <<
"," << v.y() <<
"," << v.z() <<
")";
71 double &
x,
double & y,
double &
z );
73 std::istream &
operator>>(std::istream & is, Hep3Vector & v) {
80 const Hep3Vector
HepXHat(1.0, 0.0, 0.0);
81 const Hep3Vector
HepYHat(0.0, 1.0, 0.0);
82 const Hep3Vector
HepZHat(0.0, 0.0, 1.0);
90 Hep3Vector & Hep3Vector::rotateX (
double phi1) {
91 double sinphi = std::sin(phi1);
92 double cosphi = std::cos(phi1);
94 ty = dy * cosphi - dz * sinphi;
95 dz = dz * cosphi + dy * sinphi;
100 Hep3Vector & Hep3Vector::rotateY (
double phi1) {
101 double sinphi = std::sin(phi1);
102 double cosphi = std::cos(phi1);
104 tz = dz * cosphi - dx * sinphi;
105 dx = dx * cosphi + dz * sinphi;
110 Hep3Vector & Hep3Vector::rotateZ (
double phi1) {
111 double sinphi = std::sin(phi1);
112 double cosphi = std::cos(phi1);
114 tx = dx * cosphi - dy * sinphi;
115 dy = dy * cosphi + dx * sinphi;
120 bool Hep3Vector::isNear(
const Hep3Vector & v,
double epsilon)
const {
121 double limit = dot(v)*epsilon*
epsilon;
122 return ( (*
this - v).mag2() <= limit );
125 double Hep3Vector::howNear(
const Hep3Vector & v )
const {
127 double d = (*
this - v).mag2();
129 if ( (vdv > 0) && (d < vdv) ) {
130 return std::sqrt (d/vdv);
131 }
else if ( (vdv == 0) && (d == 0) ) {
138 double Hep3Vector::deltaPhi (
const Hep3Vector & v2)
const {
139 double dphi = v2.getPhi() - getPhi();
148 double Hep3Vector::deltaR (
const Hep3Vector & v )
const {
149 double a = eta() - v.eta();
150 double b = deltaPhi(v);
151 return std::sqrt ( a*a + b*b );
154 double Hep3Vector::cosTheta(
const Hep3Vector & q)
const {
156 double ptot2 = mag2()*q.mag2();
160 arg = dot(q)/std::sqrt(ptot2);
161 if(arg > 1.0) arg = 1.0;
162 if(arg < -1.0) arg = -1.0;
167 double Hep3Vector::cos2Theta(
const Hep3Vector & q)
const {
169 double ptot2 = mag2();
170 double qtot2 = q.mag2();
171 if ( ptot2 == 0 || qtot2 == 0 ) {
175 arg = (pdq/ptot2) * (pdq/qtot2);
178 if(arg > 1.0) arg = 1.0;
183 void Hep3Vector::setEta (
double eta1) {
186 if ( (dx == 0) && (dy == 0) ) {
188 std::cerr <<
"Hep3Vector::setEta() - "
189 <<
"Attempt to set eta of zero vector -- vector is unchanged"
193 std::cerr <<
"Hep3Vector::setEta() - "
194 <<
"Attempt to set eta of vector along Z axis -- will use phi = 0"
201 double tanHalfTheta = std::exp ( -eta1 );
203 (1 - tanHalfTheta*tanHalfTheta) / (1 + tanHalfTheta*tanHalfTheta);
205 double rho1 = r1*std::sqrt(1 - cosTheta1*cosTheta1);
206 dy = rho1 * std::sin (phi1);
207 dx = rho1 * std::cos (phi1);
211 void Hep3Vector::setCylTheta (
double theta1) {
215 if ( (dx == 0) && (dy == 0) ) {
217 std::cerr <<
"Hep3Vector::setCylTheta() - "
218 <<
"Attempt to set cylTheta of zero vector -- vector is unchanged"
230 std::cerr <<
"Hep3Vector::setCylTheta() - "
231 <<
"Attempt set cylindrical theta of vector along Z axis "
232 <<
"to a non-trivial value, while keeping rho fixed -- "
233 <<
"will return zero vector" << std::endl;
237 if ( (theta1 < 0) || (theta1 >
CLHEP::pi) ) {
238 std::cerr <<
"Hep3Vector::setCylTheta() - "
239 <<
"Setting Cyl theta of a vector based on a value not in [0, PI]"
243 double phi1 (getPhi());
244 double rho1 = getRho();
245 if ( (theta1 == 0) || (theta1 ==
CLHEP::pi) ) {
246 std::cerr <<
"Hep3Vector::setCylTheta() - "
247 <<
"Attempt to set cylindrical theta to 0 or PI "
248 <<
"while keeping rho fixed -- infinite Z will be computed"
250 dz = (theta1==0) ? 1.0E72 : -1.0E72;
253 dz = rho1 / std::tan (theta1);
254 dy = rho1 * std::sin (phi1);
255 dx = rho1 * std::cos (phi1);
259 void Hep3Vector::setCylEta (
double eta1) {
263 double theta1 = 2 * std::atan ( std::exp (-eta1) );
270 if ( (dx == 0) && (dy == 0) ) {
272 std::cerr <<
"Hep3Vector::setCylEta() - "
273 <<
"Attempt to set cylEta of zero vector -- vector is unchanged"
285 std::cerr <<
"Hep3Vector::setCylEta() - "
286 <<
"Attempt set cylindrical eta of vector along Z axis "
287 <<
"to a non-trivial value, while keeping rho fixed -- "
288 <<
"will return zero vector" << std::endl;
292 double phi1 (getPhi());
293 double rho1 = getRho();
294 dz = rho1 / std::tan (theta1);
295 dy = rho1 * std::sin (phi1);
296 dx = rho1 * std::cos (phi1);
301 Hep3Vector
operator/ (
const Hep3Vector & v1,
double c ) {
307 double oneOverC = 1.0/c;
308 return Hep3Vector ( v1.x() * oneOverC,
313 Hep3Vector & Hep3Vector::operator/= (
double c) {
320 double oneOverC = 1.0/c;
HepLorentzVector operator/(const HepLorentzVector &w, double c)
const Hep3Vector HepYHat(0.0, 1.0, 0.0)
static const G4float tolerance
const Hep3Vector HepXHat(1.0, 0.0, 0.0)
std::istream & operator>>(std::istream &is, HepAxisAngle &aa)
static const double twopi
static const G4double factor
std::ostream & operator<<(std::ostream &os, const HepAxisAngle &aa)
const G4double x[NPOINTSGL]
const Hep3Vector HepZHat(0.0, 0.0, 1.0)
void ZMinput3doubles(std::istream &is, const char *type, double &x, double &y, double &z)
double epsilon(double density, double temperature)