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