| Geant4 Cross Reference |
1 // -*- C++ -*- 1 // -*- C++ -*-
>> 2 // $Id:$
2 // ------------------------------------------- 3 // ---------------------------------------------------------------------------
3 // 4 //
4 // This file is a part of the CLHEP - a Class 5 // This file is a part of the CLHEP - a Class Library for High Energy Physics.
5 // 6 //
6 // This is the implementation of that portion 7 // This is the implementation of that portion of the HepLorentzVector class
7 // which was in the original CLHEP and which d 8 // which was in the original CLHEP and which does not force loading of either
8 // Rotation.cc or LorentzRotation.cc 9 // Rotation.cc or LorentzRotation.cc
9 // 10 //
10 11
>> 12 #ifdef GNUPRAGMA
>> 13 #pragma implementation
>> 14 #endif
>> 15
11 #include "CLHEP/Vector/LorentzVector.h" 16 #include "CLHEP/Vector/LorentzVector.h"
12 17
13 #include <cmath> <<
14 #include <iostream> 18 #include <iostream>
15 19
16 namespace CLHEP { 20 namespace CLHEP {
17 21
18 double HepLorentzVector::tolerance = 22 double HepLorentzVector::tolerance =
19 Hep3Vector::ToleranceTicks * 2.22045e- 23 Hep3Vector::ToleranceTicks * 2.22045e-16;
20 double HepLorentzVector::metric = 1.0; 24 double HepLorentzVector::metric = 1.0;
21 25
22 double HepLorentzVector::operator () (int i) c 26 double HepLorentzVector::operator () (int i) const {
23 switch(i) { 27 switch(i) {
24 case X: 28 case X:
25 case Y: 29 case Y:
26 case Z: 30 case Z:
27 return pp(i); 31 return pp(i);
28 case T: 32 case T:
29 return e(); 33 return e();
30 default: 34 default:
31 std::cerr << "HepLorentzVector subscriptin 35 std::cerr << "HepLorentzVector subscripting: bad index (" << i << ")"
32 << std::endl; 36 << std::endl;
33 } 37 }
34 return 0.; 38 return 0.;
35 } 39 }
36 40
37 double & HepLorentzVector::operator () (int i) 41 double & HepLorentzVector::operator () (int i) {
38 static double dummy; 42 static double dummy;
39 switch(i) { 43 switch(i) {
40 case X: 44 case X:
41 case Y: 45 case Y:
42 case Z: 46 case Z:
43 return pp(i); 47 return pp(i);
44 case T: 48 case T:
45 return ee; 49 return ee;
46 default: 50 default:
47 std::cerr 51 std::cerr
48 << "HepLorentzVector subscripting: bad i 52 << "HepLorentzVector subscripting: bad index (" << i << ")"
49 << std::endl; 53 << std::endl;
50 return dummy; 54 return dummy;
51 } 55 }
52 } 56 }
53 57
54 HepLorentzVector & HepLorentzVector::boost 58 HepLorentzVector & HepLorentzVector::boost
55 (double bx, double by, double bz){ 59 (double bx, double by, double bz){
56 double b2 = bx*bx + by*by + bz*bz; 60 double b2 = bx*bx + by*by + bz*bz;
57 double ggamma = 1.0 / std::sqrt(1.0 - b2); 61 double ggamma = 1.0 / std::sqrt(1.0 - b2);
58 double bp = bx*x() + by*y() + bz*z(); 62 double bp = bx*x() + by*y() + bz*z();
59 double gamma2 = b2 > 0 ? (ggamma - 1.0)/b2 : 63 double gamma2 = b2 > 0 ? (ggamma - 1.0)/b2 : 0.0;
60 64
61 setX(x() + gamma2*bp*bx + ggamma*bx*t()); 65 setX(x() + gamma2*bp*bx + ggamma*bx*t());
62 setY(y() + gamma2*bp*by + ggamma*by*t()); 66 setY(y() + gamma2*bp*by + ggamma*by*t());
63 setZ(z() + gamma2*bp*bz + ggamma*bz*t()); 67 setZ(z() + gamma2*bp*bz + ggamma*bz*t());
64 setT(ggamma*(t() + bp)); 68 setT(ggamma*(t() + bp));
65 return *this; 69 return *this;
66 } 70 }
67 71
68 HepLorentzVector & HepLorentzVector::rotateX(d 72 HepLorentzVector & HepLorentzVector::rotateX(double a) {
69 pp.rotateX(a); 73 pp.rotateX(a);
70 return *this; 74 return *this;
71 } 75 }
72 HepLorentzVector & HepLorentzVector::rotateY(d 76 HepLorentzVector & HepLorentzVector::rotateY(double a) {
73 pp.rotateY(a); 77 pp.rotateY(a);
74 return *this; 78 return *this;
75 } 79 }
76 HepLorentzVector & HepLorentzVector::rotateZ(d 80 HepLorentzVector & HepLorentzVector::rotateZ(double a) {
77 pp.rotateZ(a); 81 pp.rotateZ(a);
78 return *this; 82 return *this;
79 } 83 }
80 84
81 HepLorentzVector & HepLorentzVector::rotateUz( 85 HepLorentzVector & HepLorentzVector::rotateUz(const Hep3Vector &v1) {
82 pp.rotateUz(v1); 86 pp.rotateUz(v1);
83 return *this; 87 return *this;
84 } 88 }
85 89
86 std::ostream & operator<< (std::ostream & os, 90 std::ostream & operator<< (std::ostream & os, const HepLorentzVector & v1)
87 { 91 {
88 return os << "(" << v1.x() << "," << v1.y() 92 return os << "(" << v1.x() << "," << v1.y() << "," << v1.z()
89 << ";" << v1.t() << ")"; 93 << ";" << v1.t() << ")";
90 } 94 }
91 95
92 std::istream & operator>> (std::istream & is, 96 std::istream & operator>> (std::istream & is, HepLorentzVector & v1) {
93 97
94 // Required format is ( a, b, c; d ) that is, 98 // Required format is ( a, b, c; d ) that is, four numbers, preceded by
95 // (, followed by ), components of the spatial 99 // (, followed by ), components of the spatial vector separated by commas,
96 // time component separated by semicolon. The 100 // time component separated by semicolon. The four numbers are taken
97 // as x, y, z, t. 101 // as x, y, z, t.
98 102
99 double x, y, z, t; 103 double x, y, z, t;
100 char c; 104 char c;
101 105
102 is >> std::ws >> c; 106 is >> std::ws >> c;
103 // ws is defined to invoke eatwhite(istrea 107 // ws is defined to invoke eatwhite(istream & )
104 // see (Stroustrup gray book) page 333 and 108 // see (Stroustrup gray book) page 333 and 345.
105 if (is.fail() || c != '(' ) { 109 if (is.fail() || c != '(' ) {
106 std::cerr << "Could not find required open 110 std::cerr << "Could not find required opening parenthesis "
107 << "in input of a HepLorentzVector" << 111 << "in input of a HepLorentzVector" << std::endl;
108 return is; 112 return is;
109 } 113 }
110 114
111 is >> x >> std::ws >> c; 115 is >> x >> std::ws >> c;
112 if (is.fail() || c != ',' ) { 116 if (is.fail() || c != ',' ) {
113 std::cerr << "Could not find x value and r 117 std::cerr << "Could not find x value and required trailing comma "
114 << "in input of a HepLorentzVector" << 118 << "in input of a HepLorentzVector" << std::endl;
115 return is; 119 return is;
116 } 120 }
117 121
118 is >> y >> std::ws >> c; 122 is >> y >> std::ws >> c;
119 if (is.fail() || c != ',' ) { 123 if (is.fail() || c != ',' ) {
120 std::cerr << "Could not find y value and r 124 std::cerr << "Could not find y value and required trailing comma "
121 << "in input of a HepLorentzVec 125 << "in input of a HepLorentzVector" << std::endl;
122 return is; 126 return is;
123 } 127 }
124 128
125 is >> z >> std::ws >> c; 129 is >> z >> std::ws >> c;
126 if (is.fail() || c != ';' ) { 130 if (is.fail() || c != ';' ) {
127 std::cerr << "Could not find z value and r 131 std::cerr << "Could not find z value and required trailing semicolon "
128 << "in input of a HepLorentzVector" << s 132 << "in input of a HepLorentzVector" << std::endl;
129 return is; 133 return is;
130 } 134 }
131 135
132 is >> t >> std::ws >> c; 136 is >> t >> std::ws >> c;
133 if (is.fail() || c != ')' ) { 137 if (is.fail() || c != ')' ) {
134 std::cerr << "Could not find t value and r 138 std::cerr << "Could not find t value and required close parenthesis "
135 << "in input of a HepLorentzVector" << st 139 << "in input of a HepLorentzVector" << std::endl;
136 return is; 140 return is;
137 } 141 }
138 142
139 v1.setX(x); 143 v1.setX(x);
140 v1.setY(y); 144 v1.setY(y);
141 v1.setZ(z); 145 v1.setZ(z);
142 v1.setT(t); 146 v1.setT(t);
143 return is; 147 return is;
144 } 148 }
145 149
146 // The following were added when ZOOM classes 150 // The following were added when ZOOM classes were merged in:
147 151
148 HepLorentzVector & HepLorentzVector::operator 152 HepLorentzVector & HepLorentzVector::operator /= (double c) {
149 // if (c == 0) { 153 // if (c == 0) {
150 // std::cerr << "HepLorentzVector::operator 154 // std::cerr << "HepLorentzVector::operator /=() - "
151 // << "Attempt to do LorentzVector /= 0 - 155 // << "Attempt to do LorentzVector /= 0 -- \n"
152 // << "division by zero would produce inf 156 // << "division by zero would produce infinite or NAN components"
153 // << std::endl; 157 // << std::endl;
154 // } 158 // }
155 double oneOverC = 1.0/c; 159 double oneOverC = 1.0/c;
156 pp *= oneOverC; 160 pp *= oneOverC;
157 ee *= oneOverC; 161 ee *= oneOverC;
158 return *this; 162 return *this;
159 } /* w /= c */ 163 } /* w /= c */
160 164
161 HepLorentzVector operator / (const HepLorentzV 165 HepLorentzVector operator / (const HepLorentzVector & w, double c) {
162 // if (c == 0) { 166 // if (c == 0) {
163 // std::cerr << "HepLorentzVector::operator 167 // std::cerr << "HepLorentzVector::operator /() - "
164 // << "Attempt to do LorentzVector / 0 -- 168 // << "Attempt to do LorentzVector / 0 -- \n"
165 // << "division by zero would produce inf 169 // << "division by zero would produce infinite or NAN components"
166 // << std::endl; 170 // << std::endl;
167 // } 171 // }
168 double oneOverC = 1.0/c; 172 double oneOverC = 1.0/c;
169 return HepLorentzVector (w.getV() * oneOverC 173 return HepLorentzVector (w.getV() * oneOverC,
170 w.getT() * oneOverC); 174 w.getT() * oneOverC);
171 } /* LV = w / c */ 175 } /* LV = w / c */
172 176
173 Hep3Vector HepLorentzVector::boostVector() con 177 Hep3Vector HepLorentzVector::boostVector() const {
174 if (ee == 0) { 178 if (ee == 0) {
175 if (pp.mag2() == 0) { 179 if (pp.mag2() == 0) {
176 return Hep3Vector(0,0,0); 180 return Hep3Vector(0,0,0);
177 } else { 181 } else {
178 std::cerr << "HepLorentzVector::boostVec 182 std::cerr << "HepLorentzVector::boostVector() - "
179 << "boostVector computed for LorentzVe 183 << "boostVector computed for LorentzVector with t=0 -- infinite result"
180 << std::endl; 184 << std::endl;
181 return pp/ee; 185 return pp/ee;
182 } 186 }
183 } 187 }
184 if (restMass2() <= 0) { 188 if (restMass2() <= 0) {
185 std::cerr << "HepLorentzVector::boostVecto 189 std::cerr << "HepLorentzVector::boostVector() - "
186 << "boostVector computed for a non-timel 190 << "boostVector computed for a non-timelike LorentzVector " << std::endl;
187 // result will make analytic sense but 191 // result will make analytic sense but is physically meaningless
188 } 192 }
189 return pp * (1./ee); 193 return pp * (1./ee);
190 } /* boostVector */ 194 } /* boostVector */
191 195
192 196
193 HepLorentzVector & HepLorentzVector::boostX (d 197 HepLorentzVector & HepLorentzVector::boostX (double bbeta){
194 double b2 = bbeta*bbeta; 198 double b2 = bbeta*bbeta;
195 if (b2 >= 1) { 199 if (b2 >= 1) {
196 std::cerr << "HepLorentzVector::boostX() - 200 std::cerr << "HepLorentzVector::boostX() - "
197 << "boost along X with beta >= 1 (speed 201 << "boost along X with beta >= 1 (speed of light) -- \n"
198 << "no boost done" << std::endl; 202 << "no boost done" << std::endl;
199 } else { 203 } else {
200 double ggamma = std::sqrt(1./(1-b2)); 204 double ggamma = std::sqrt(1./(1-b2));
201 double tt = ee; 205 double tt = ee;
202 ee = ggamma*(ee + bbeta*pp.getX()); 206 ee = ggamma*(ee + bbeta*pp.getX());
203 pp.setX(ggamma*(pp.getX() + bbeta*tt)); 207 pp.setX(ggamma*(pp.getX() + bbeta*tt));
204 } 208 }
205 return *this; 209 return *this;
206 } /* boostX */ 210 } /* boostX */
207 211
208 HepLorentzVector & HepLorentzVector::boostY (d 212 HepLorentzVector & HepLorentzVector::boostY (double bbeta){
209 double b2 = bbeta*bbeta; 213 double b2 = bbeta*bbeta;
210 if (b2 >= 1) { 214 if (b2 >= 1) {
211 std::cerr << "HepLorentzVector::boostY() - 215 std::cerr << "HepLorentzVector::boostY() - "
212 << "boost along Y with beta >= 1 (speed 216 << "boost along Y with beta >= 1 (speed of light) -- \n"
213 << "no boost done" << std::endl; 217 << "no boost done" << std::endl;
214 } else { 218 } else {
215 double ggamma = std::sqrt(1./(1-b2)); 219 double ggamma = std::sqrt(1./(1-b2));
216 double tt = ee; 220 double tt = ee;
217 ee = ggamma*(ee + bbeta*pp.getY()); 221 ee = ggamma*(ee + bbeta*pp.getY());
218 pp.setY(ggamma*(pp.getY() + bbeta*tt)); 222 pp.setY(ggamma*(pp.getY() + bbeta*tt));
219 } 223 }
220 return *this; 224 return *this;
221 } /* boostY */ 225 } /* boostY */
222 226
223 HepLorentzVector & HepLorentzVector::boostZ (d 227 HepLorentzVector & HepLorentzVector::boostZ (double bbeta){
224 double b2 = bbeta*bbeta; 228 double b2 = bbeta*bbeta;
225 if (b2 >= 1) { 229 if (b2 >= 1) {
226 std::cerr << "HepLorentzVector::boostZ() - 230 std::cerr << "HepLorentzVector::boostZ() - "
227 << "boost along Z with beta >= 1 (speed 231 << "boost along Z with beta >= 1 (speed of light) -- \n"
228 << "no boost done" << std::endl; 232 << "no boost done" << std::endl;
229 } else { 233 } else {
230 double ggamma = std::sqrt(1./(1-b2)); 234 double ggamma = std::sqrt(1./(1-b2));
231 double tt = ee; 235 double tt = ee;
232 ee = ggamma*(ee + bbeta*pp.getZ()); 236 ee = ggamma*(ee + bbeta*pp.getZ());
233 pp.setZ(ggamma*(pp.getZ() + bbeta*tt)); 237 pp.setZ(ggamma*(pp.getZ() + bbeta*tt));
234 } 238 }
235 return *this; 239 return *this;
236 } /* boostZ */ 240 } /* boostZ */
237 241
238 double HepLorentzVector::setTolerance ( double 242 double HepLorentzVector::setTolerance ( double tol ) {
239 // Set the tolerance for two LorentzVectors to 243 // Set the tolerance for two LorentzVectors to be considered near each other
240 double oldTolerance (tolerance); 244 double oldTolerance (tolerance);
241 tolerance = tol; 245 tolerance = tol;
242 return oldTolerance; 246 return oldTolerance;
243 } 247 }
244 248
245 double HepLorentzVector::getTolerance ( ) { 249 double HepLorentzVector::getTolerance ( ) {
246 // Get the tolerance for two LorentzVectors to 250 // Get the tolerance for two LorentzVectors to be considered near each other
247 return tolerance; 251 return tolerance;
248 } 252 }
249 253
250 } // namespace CLHEP 254 } // namespace CLHEP
251 255