CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

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