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