| Geant4 Cross Reference |
1 // -*- C++ -*- 1
2 // -------------------------------------------
3 //
4 // This file is a part of the CLHEP - a Class
5 //
6 // This is the implementation of that portion
7 // which was in the original CLHEP and which d
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-
20 double HepLorentzVector::metric = 1.0;
21
22 double HepLorentzVector::operator () (int i) c
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 subscriptin
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 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 :
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(d
69 pp.rotateX(a);
70 return *this;
71 }
72 HepLorentzVector & HepLorentzVector::rotateY(d
73 pp.rotateY(a);
74 return *this;
75 }
76 HepLorentzVector & HepLorentzVector::rotateZ(d
77 pp.rotateZ(a);
78 return *this;
79 }
80
81 HepLorentzVector & HepLorentzVector::rotateUz(
82 pp.rotateUz(v1);
83 return *this;
84 }
85
86 std::ostream & operator<< (std::ostream & os,
87 {
88 return os << "(" << v1.x() << "," << v1.y()
89 << ";" << v1.t() << ")";
90 }
91
92 std::istream & operator>> (std::istream & is,
93
94 // Required format is ( a, b, c; d ) that is,
95 // (, followed by ), components of the spatial
96 // time component separated by semicolon. The
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(istrea
104 // see (Stroustrup gray book) page 333 and
105 if (is.fail() || c != '(' ) {
106 std::cerr << "Could not find required open
107 << "in input of a HepLorentzVector" <<
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 r
114 << "in input of a HepLorentzVector" <<
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 r
121 << "in input of a HepLorentzVec
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 r
128 << "in input of a HepLorentzVector" << s
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 r
135 << "in input of a HepLorentzVector" << st
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
147
148 HepLorentzVector & HepLorentzVector::operator
149 // if (c == 0) {
150 // std::cerr << "HepLorentzVector::operator
151 // << "Attempt to do LorentzVector /= 0 -
152 // << "division by zero would produce inf
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 HepLorentzV
162 // if (c == 0) {
163 // std::cerr << "HepLorentzVector::operator
164 // << "Attempt to do LorentzVector / 0 --
165 // << "division by zero would produce inf
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() con
174 if (ee == 0) {
175 if (pp.mag2() == 0) {
176 return Hep3Vector(0,0,0);
177 } else {
178 std::cerr << "HepLorentzVector::boostVec
179 << "boostVector computed for LorentzVe
180 << std::endl;
181 return pp/ee;
182 }
183 }
184 if (restMass2() <= 0) {
185 std::cerr << "HepLorentzVector::boostVecto
186 << "boostVector computed for a non-timel
187 // result will make analytic sense but
188 }
189 return pp * (1./ee);
190 } /* boostVector */
191
192
193 HepLorentzVector & HepLorentzVector::boostX (d
194 double b2 = bbeta*bbeta;
195 if (b2 >= 1) {
196 std::cerr << "HepLorentzVector::boostX() -
197 << "boost along X with beta >= 1 (speed
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 (d
209 double b2 = bbeta*bbeta;
210 if (b2 >= 1) {
211 std::cerr << "HepLorentzVector::boostY() -
212 << "boost along Y with beta >= 1 (speed
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 (d
224 double b2 = bbeta*bbeta;
225 if (b2 >= 1) {
226 std::cerr << "HepLorentzVector::boostZ() -
227 << "boost along Z with beta >= 1 (speed
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
239 // Set the tolerance for two LorentzVectors to
240 double oldTolerance (tolerance);
241 tolerance = tol;
242 return oldTolerance;
243 }
244
245 double HepLorentzVector::getTolerance ( ) {
246 // Get the tolerance for two LorentzVectors to
247 return tolerance;
248 }
249
250 } // namespace CLHEP
251