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