| 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 // SpaceVector 6 // SpaceVector
7 // 7 //
8 // This is the implementation of the subset of 8 // This is the implementation of the subset of those methods of the Hep3Vector
9 // class which originated from the ZOOM SpaceV 9 // class which originated from the ZOOM SpaceVector class *and* which involve
10 // intrinsic properties or propeties relative 10 // intrinsic properties or propeties relative to a second vector.
11 // 11 //
12 12
>> 13 #ifdef GNUPRAGMA
>> 14 #pragma implementation
>> 15 #endif
>> 16
13 #include "CLHEP/Vector/ThreeVector.h" 17 #include "CLHEP/Vector/ThreeVector.h"
14 18
15 #include <cmath> 19 #include <cmath>
16 20
17 namespace CLHEP { 21 namespace CLHEP {
18 22
19 //-******************************** 23 //-********************************
20 // - 5 - 24 // - 5 -
21 // Intrinsic properties of a vector 25 // Intrinsic properties of a vector
22 // and properties relative to a direction 26 // and properties relative to a direction
23 // 27 //
24 //-******************************** 28 //-********************************
25 29
26 double Hep3Vector::beta() const { 30 double Hep3Vector::beta() const {
27 double b = std::sqrt(mag2()); 31 double b = std::sqrt(mag2());
28 // if (b >= 1) { 32 // if (b >= 1) {
29 // std::cerr << "Hep3Vector::beta() - " 33 // std::cerr << "Hep3Vector::beta() - "
30 // << "Beta taken for Hep3Vector of at le 34 // << "Beta taken for Hep3Vector of at least unit length" << std::endl;
31 // } 35 // }
32 return b; 36 return b;
33 } 37 }
34 38
35 double Hep3Vector::gamma() const { 39 double Hep3Vector::gamma() const {
36 double bbeta = std::sqrt(mag2()); 40 double bbeta = std::sqrt(mag2());
37 // if (bbeta == 1) { 41 // if (bbeta == 1) {
38 // std::cerr << "Hep3Vector::gamma() - " 42 // std::cerr << "Hep3Vector::gamma() - "
39 // << "Gamma taken for Hep3Vector of unit 43 // << "Gamma taken for Hep3Vector of unit magnitude -- infinite result"
40 // << std::endl; 44 // << std::endl;
41 // } 45 // }
42 // if (bbeta > 1) { 46 // if (bbeta > 1) {
43 // std::cerr << "Hep3Vector::gamma() - " 47 // std::cerr << "Hep3Vector::gamma() - "
44 // << "Gamma taken for Hep3Vector of more 48 // << "Gamma taken for Hep3Vector of more than unit magnitude -- \n"
45 // << "the sqrt function would return NAN 49 // << "the sqrt function would return NAN" << std::endl;
46 // } 50 // }
47 return 1/std::sqrt(1-bbeta*bbeta); 51 return 1/std::sqrt(1-bbeta*bbeta);
48 } 52 }
49 53
50 double Hep3Vector::rapidity() const { 54 double Hep3Vector::rapidity() const {
51 // if (std::fabs(z()) == 1) { << 55 // if (std::fabs(dz) == 1) {
52 // std::cerr << "Hep3Vector::rapidity() - " 56 // std::cerr << "Hep3Vector::rapidity() - "
53 // << "Rapidity in Z direction taken for 57 // << "Rapidity in Z direction taken for Hep3Vector with |Z| = 1 -- \n"
54 // << "the log should return infinity" <, 58 // << "the log should return infinity" <, std::endl;
55 // } 59 // }
56 // if (std::fabs(z()) > 1) { << 60 // if (std::fabs(dz) > 1) {
57 // std::cerr << "Hep3Vector::rapidity() - " 61 // std::cerr << "Hep3Vector::rapidity() - "
58 // << "Rapidity in Z direction taken for 62 // << "Rapidity in Z direction taken for Hep3Vector with |Z| > 1 -- \n"
59 // << "the log would return a NAN" << std 63 // << "the log would return a NAN" << std::endl;
60 // } 64 // }
61 // Want inverse std::tanh(dz): 65 // Want inverse std::tanh(dz):
62 return (.5 * std::log((1+z())/(1-z())) ); << 66 return (.5 * std::log((1+dz)/(1-dz)) );
63 } 67 }
64 68
65 double Hep3Vector::coLinearRapidity() const { 69 double Hep3Vector::coLinearRapidity() const {
66 double b = beta(); 70 double b = beta();
67 // if (b == 1) { 71 // if (b == 1) {
68 // std::cerr << "Hep3Vector::coLinearRapidi 72 // std::cerr << "Hep3Vector::coLinearRapidity() - "
69 // << "Co-linear Rapidity taken for Hep3V 73 // << "Co-linear Rapidity taken for Hep3Vector of unit length -- \n"
70 // << "the log should return infinity" << 74 // << "the log should return infinity" << std::endl;
71 // } 75 // }
72 // if (b > 1) { 76 // if (b > 1) {
73 // std::cerr << "Hep3Vector::coLinearRapidi 77 // std::cerr << "Hep3Vector::coLinearRapidity() - "
74 // << "Co-linear Rapidity taken for Hep3V 78 // << "Co-linear Rapidity taken for Hep3Vector of more than unit length -- \n"
75 // << "the log would return a NAN" << std 79 // << "the log would return a NAN" << std::endl;
76 // } 80 // }
77 // Want inverse std::tanh(b): 81 // Want inverse std::tanh(b):
78 return (.5 * std::log((1+b)/(1-b)) ); 82 return (.5 * std::log((1+b)/(1-b)) );
79 } 83 }
80 84
81 //-******************************************* 85 //-***********************************************
82 // Other properties relative to a reference ve 86 // Other properties relative to a reference vector
83 //-******************************************* 87 //-***********************************************
84 88
85 Hep3Vector Hep3Vector::project (const Hep3Vect 89 Hep3Vector Hep3Vector::project (const Hep3Vector & v2) const {
86 double mag2v2 = v2.mag2(); 90 double mag2v2 = v2.mag2();
87 if (mag2v2 == 0) { 91 if (mag2v2 == 0) {
88 std::cerr << "Hep3Vector::project() - " 92 std::cerr << "Hep3Vector::project() - "
89 << "Attempt to take projection of vector 93 << "Attempt to take projection of vector against zero reference vector"
90 << std::endl; 94 << std::endl;
91 return project(); 95 return project();
92 } 96 }
93 return ( v2 * (dot(v2)/mag2v2) ); 97 return ( v2 * (dot(v2)/mag2v2) );
94 } 98 }
95 99
96 double Hep3Vector::rapidity(const Hep3Vector & 100 double Hep3Vector::rapidity(const Hep3Vector & v2) const {
97 double vmag = v2.mag(); 101 double vmag = v2.mag();
98 if ( vmag == 0 ) { 102 if ( vmag == 0 ) {
99 std::cerr << "Hep3Vector::rapidity() - " 103 std::cerr << "Hep3Vector::rapidity() - "
100 << "Rapidity taken with respect to zero 104 << "Rapidity taken with respect to zero vector" << std::endl;
101 return 0; 105 return 0;
102 } 106 }
103 double z1 = dot(v2)/vmag; 107 double z1 = dot(v2)/vmag;
104 // if (std::fabs(z1) >= 1) { 108 // if (std::fabs(z1) >= 1) {
105 // std::cerr << "Hep3Vector::rapidity() - " 109 // std::cerr << "Hep3Vector::rapidity() - "
106 // << "Rapidity taken for too large a Hep 110 // << "Rapidity taken for too large a Hep3Vector "
107 // << "-- would return infinity or NAN" < 111 // << "-- would return infinity or NAN" << std::endl;
108 // } 112 // }
109 // Want inverse std::tanh(z): 113 // Want inverse std::tanh(z):
110 return (.5 * std::log((1+z1)/(1-z1)) ); 114 return (.5 * std::log((1+z1)/(1-z1)) );
111 } 115 }
112 116
113 double Hep3Vector::eta(const Hep3Vector & v2) 117 double Hep3Vector::eta(const Hep3Vector & v2) const {
114 // Defined as -std::log ( std::tan ( .5* 118 // Defined as -std::log ( std::tan ( .5* theta(u) ) );
115 // 119 //
116 // Quicker is to use cosTheta: 120 // Quicker is to use cosTheta:
117 // std::tan (theta/2) = std::sin(theta)/(1 + 121 // std::tan (theta/2) = std::sin(theta)/(1 + std::cos(theta))
118 122
119 double r1 = getR(); 123 double r1 = getR();
120 double v2r = v2.mag(); 124 double v2r = v2.mag();
121 if ( (r1 == 0) || (v2r == 0) ) { 125 if ( (r1 == 0) || (v2r == 0) ) {
122 std::cerr << "Hep3Vector::eta() - " 126 std::cerr << "Hep3Vector::eta() - "
123 << "Cannot find pseudorapidity of a zero 127 << "Cannot find pseudorapidity of a zero vector relative to a vector"
124 << std::endl; 128 << std::endl;
125 return 0.; 129 return 0.;
126 } 130 }
127 double c = dot(v2)/(r1*v2r); 131 double c = dot(v2)/(r1*v2r);
128 if ( c >= 1 ) { 132 if ( c >= 1 ) {
129 c = 1; //-| We don't want to return NAN b 133 c = 1; //-| We don't want to return NAN because of roundoff
130 std::cerr << "Hep3Vector::eta() - " 134 std::cerr << "Hep3Vector::eta() - "
131 << "Pseudorapidity of vector relative to 135 << "Pseudorapidity of vector relative to parallel vector -- \n"
132 << "will give infinite result" << std::e 136 << "will give infinite result" << std::endl;
133 // We can just go on; tangent will 137 // We can just go on; tangent will be 0, so
134 // std::log (tangent) will be -INFIN 138 // std::log (tangent) will be -INFINITY, so result
135 // will be +INFINITY. 139 // will be +INFINITY.
136 } 140 }
137 if ( c <= -1 ) { 141 if ( c <= -1 ) {
138 std::cerr << "Hep3Vector::eta() - " 142 std::cerr << "Hep3Vector::eta() - "
139 << "Pseudorapidity of vector relative to 143 << "Pseudorapidity of vector relative to anti-parallel vector -- \n"
140 << "will give negative infinite result"< 144 << "will give negative infinite result"<< std::endl;
141 //-| We don't want to return NAN bec 145 //-| We don't want to return NAN because of roundoff
142 return ( negativeInfinity() ); 146 return ( negativeInfinity() );
143 // If we just went on, the tangent 147 // If we just went on, the tangent would be NAN
144 // so return would be NAN. But the 148 // so return would be NAN. But the proper limit
145 // of tan is +Infinity, so the retur 149 // of tan is +Infinity, so the return should be
146 // -INFINITY. 150 // -INFINITY.
147 } 151 }
148 152
149 double tangent = std::sqrt (1-c*c) / ( 1 + c 153 double tangent = std::sqrt (1-c*c) / ( 1 + c );
150 return (- std::log (tangent)); 154 return (- std::log (tangent));
151 155
152 } /* eta (u) */ 156 } /* eta (u) */
153 157
154 158
155 } // namespace CLHEP 159 } // namespace CLHEP
156 160