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 
22 namespace CLHEP {
23 
24 double 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 
88 std::ostream & operator<< (std::ostream & os, const HepLorentzVector & v1)
89 {
90  return os << "(" << v1.x() << "," << v1.y() << "," << v1.z()
91  << ";" << v1.t() << ")";
92 }
93 
94 std::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 
163 if (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 
234 double 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 
246 double HepLorentzVector::tolerance =
247  Hep3Vector::ToleranceTicks * 2.22045e-16;
248 double HepLorentzVector::metric = 1.0;
249 
250 
251 } // namespace CLHEP
HepLorentzVector & boostX(double beta)
double getZ() const
static double getTolerance()
HepLorentzVector operator/(const HepLorentzVector &, double a)
HepLorentzVector & boostY(double beta)
HepLorentzVector & rotateX(double)
double restMass2() const
Hep3Vector & rotateX(double)
Definition: ThreeVector.cc:124
Hep3Vector boostVector() const
double getX() const
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:72
#define ZMthrowC(A)
void setX(double)
double mag2() const
static double setTolerance(double tol)
double getY() const
void setZ(double)
HepLorentzVector & rotateUz(const Hep3Vector &)
HepLorentzVector & operator/=(double)
HepLorentzVector & rotateY(double)
double operator()(int) const
#define ZMthrowA(A)
std::istream & operator>>(std::istream &is, HepAxisAngle &aa)
Definition: AxisAngle.cc:96
Hep3Vector getV() const
Hep3Vector & rotateZ(double)
Definition: ThreeVector.cc:144
HepLorentzVector & rotateZ(double)
HepLorentzVector & boostZ(double beta)
Hep3Vector & rotateY(double)
Definition: ThreeVector.cc:134
std::ostream & operator<<(std::ostream &os, const HepAxisAngle &aa)
Definition: AxisAngle.cc:86
void setY(double)
HepLorentzVector & boost(double, double, double)