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