Geant4  10.02.p01
LorentzVectorK.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 // ---------------------------------------------------------------------------
3 //
4 // This file is a part of the CLHEP - a Class Library for High Energy Physics.
5 //
6 // This is part of the implementation of the HepLorentzVector class:
7 // Those methods which originated from ZOOM and which deal with relativistic
8 // kinematic properties.
9 //
10 
11 #ifdef GNUPRAGMA
12 #pragma implementation
13 #endif
14 
15 #include "CLHEP/Vector/LorentzVector.h"
16 
17 #include <cmath>
18 
19 namespace CLHEP {
20 
21 //-******************
22 // Metric flexibility
23 //-******************
24 
25 ZMpvMetric_t HepLorentzVector::setMetric( ZMpvMetric_t a1 ) {
26  ZMpvMetric_t oldMetric = (metric > 0) ? TimePositive : TimeNegative;
27  if ( a1 == TimeNegative ) {
28  metric = -1.0;
29  } else {
30  metric = 1.0;
31  }
32  return oldMetric;
33 }
34 
35 ZMpvMetric_t HepLorentzVector::getMetric() {
36  return ( (metric > 0) ? TimePositive : TimeNegative );
37 }
38 
39 //-********
40 // plus
41 // minus
42 //-********
43 
44 double HepLorentzVector::plus (const Hep3Vector & ref) const {
45  double r = ref.mag();
46  if (r == 0) {
47  std::cerr << "HepLorentzVector::plus() - "
48  << "A zero vector used as reference to LorentzVector plus-part"
49  << std::endl;
50  return ee;
51  }
52  return ee + pp.dot(ref)/r;
53 } /* plus */
54 
55 double HepLorentzVector::minus (const Hep3Vector & ref) const {
56  double r = ref.mag();
57  if (r == 0) {
58  std::cerr << "HepLorentzVector::minus() - "
59  << "A zero vector used as reference to LorentzVector minus-part"
60  << std::endl;
61  return ee;
62  }
63  return ee - pp.dot(ref)/r;
64 } /* plus */
65 
66 HepLorentzVector HepLorentzVector::rest4Vector() const {
67  return HepLorentzVector (0, 0, 0, (t() < 0.0 ? -m() : m()));
68 }
69 
70 //-********
71 // beta
72 // gamma
73 //-********
74 
75 double HepLorentzVector::beta() const {
76  if (ee == 0) {
77  if (pp.mag2() == 0) {
78  return 0;
79  } else {
80  std::cerr << "HepLorentzVector::beta() - "
81  << "beta computed for HepLorentzVector with t=0 -- infinite result"
82  << std::endl;
83  return 1./ee;
84  }
85  }
86 // if (restMass2() <= 0) {
87 // std::cerr << "HepLorentzVector::beta() - "
88 // << "beta computed for a non-timelike HepLorentzVector" << std::endl;
89 // // result will make analytic sense but is physically meaningless
90 // }
91  return std::sqrt (pp.mag2() / (ee*ee)) ;
92 } /* beta */
93 
94 double HepLorentzVector::gamma() const {
95  double v2 = pp.mag2();
96  double t2 = ee*ee;
97  if (ee == 0) {
98  if (pp.mag2() == 0) {
99  return 1;
100  } else {
101  std::cerr << "HepLorentzVector::gamma() - "
102  << "gamma computed for HepLorentzVector with t=0 -- zero result"
103  << std::endl;
104  return 0;
105  }
106  }
107  if (t2 < v2) {
108  std::cerr << "HepLorentzVector::gamma() - "
109  << "gamma computed for a spacelike HepLorentzVector -- imaginary result"
110  << std::endl;
111  // analytic result would be imaginary.
112  return 0;
113 // } else if ( t2 == v2 ) {
114 // std::cerr << "HepLorentzVector::gamma() - "
115 // << "gamma computed for a lightlike HepLorentzVector -- infinite result"
116 // << std::endl;
117  }
118  return 1./std::sqrt(1. - v2/t2 );
119 } /* gamma */
120 
121 
122 //-***************
123 // rapidity
124 // pseudorapidity
125 // eta
126 //-***************
127 
128 double HepLorentzVector::rapidity() const {
129  double z1 = pp.getZ();
130 // if (std::fabs(ee) == std::fabs(z1)) {
131 // std::cerr << "HepLorentzVector::rapidity() - "
132 // << "rapidity for 4-vector with |E| = |Pz| -- infinite result"
133 // << std::endl;
134 // }
135  if (std::fabs(ee) < std::fabs(z1)) {
136  std::cerr << "HepLorentzVector::rapidity() - "
137  << "rapidity for spacelike 4-vector with |E| < |Pz| -- undefined"
138  << std::endl;
139  return 0;
140  }
141  double q = (ee + z1) / (ee - z1);
142  //-| This cannot be negative now, since both numerator
143  //-| and denominator have the same sign as ee.
144  return .5 * std::log(q);
145 } /* rapidity */
146 
147 double HepLorentzVector::rapidity(const Hep3Vector & ref) const {
148  double r = ref.mag2();
149  if (r == 0) {
150  std::cerr << "HepLorentzVector::rapidity() - "
151  << "A zero vector used as reference to LorentzVector rapidity"
152  << std::endl;
153  return 0;
154  }
155  double vdotu = pp.dot(ref)/std::sqrt(r);
156 // if (std::fabs(ee) == std::fabs(vdotu)) {
157 // std::cerr << "HepLorentzVector::rapidity() - "
158 // << "rapidity for 4-vector with |E| = |Pu| -- infinite result"
159 // << std::endl;
160 // }
161  if (std::fabs(ee) < std::fabs(vdotu)) {
162  std::cerr << "HepLorentzVector::rapidity() - "
163  << "rapidity for spacelike 4-vector with |E| < |P*ref| -- undefined "
164  << std::endl;
165  return 0;
166  }
167  double q = (ee + vdotu) / (ee - vdotu);
168  return .5 * std::log(q);
169 } /* rapidity(ref) */
170 
171 double HepLorentzVector::coLinearRapidity() const {
172  double v1 = pp.mag();
173 // if (std::fabs(ee) == std::fabs(v1)) {
174 // std::cerr << "HepLorentzVector::coLinearRapidity() - "
175 // << "co-Linear rapidity for 4-vector with |E| = |P| -- infinite result"
176 // << std::endl;
177 // }
178  if (std::fabs(ee) < std::fabs(v1)) {
179  std::cerr << "HepLorentzVector::coLinearRapidity() - "
180  << "co-linear rapidity for spacelike 4-vector -- undefined"
181  << std::endl;
182  return 0;
183  }
184  double q = (ee + v1) / (ee - v1);
185  return .5 * std::log(q);
186 } /* rapidity */
187 
188 //-*************
189 // invariantMass
190 //-*************
191 
192 double HepLorentzVector::invariantMass(const HepLorentzVector & w) const {
193  double m1 = invariantMass2(w);
194  if (m1 < 0) {
195  // We should find out why:
196  if ( ee * w.ee < 0 ) {
197  std::cerr << "HepLorentzVector::invariantMass() - "
198  << "invariant mass meaningless: \n"
199  << "a negative-mass input led to spacelike 4-vector sum" << std::endl;
200  return 0;
201  } else if ( (isSpacelike() && !isLightlike()) ||
202  (w.isSpacelike() && !w.isLightlike()) ) {
203  std::cerr << "HepLorentzVector::invariantMass() - "
204  << "invariant mass meaningless because of spacelike input"
205  << std::endl;
206  return 0;
207  } else {
208  // Invariant mass squared for a pair of timelike or lightlike vectors
209  // mathematically cannot be negative. If the vectors are within the
210  // tolerance of being lightlike or timelike, we can assume that prior
211  // or current roundoffs have caused the negative result, and return 0
212  // without comment.
213  return 0;
214  }
215  }
216  return (ee+w.ee >=0 ) ? std::sqrt(m1) : - std::sqrt(m1);
217 } /* invariantMass */
218 
219 //-***************
220 // findBoostToCM
221 //-***************
222 
223 Hep3Vector HepLorentzVector::findBoostToCM() const {
224  return -boostVector();
225 } /* boostToCM() */
226 
227 Hep3Vector HepLorentzVector::findBoostToCM (const HepLorentzVector & w) const {
228  double t1 = ee + w.ee;
229  Hep3Vector v1 = pp + w.pp;
230  if (t1 == 0) {
231  if (v1.mag2() == 0) {
232  return Hep3Vector(0,0,0);
233  } else {
234  std::cerr << "HepLorentzVector::findBoostToCM() - "
235  << "boostToCM computed for two 4-vectors with combined t=0 -- "
236  << "infinite result" << std::endl;
237  return Hep3Vector(v1*(1./t1)); // Yup, 1/0 -- that is how we return infinity
238  }
239  }
240 // if (t1*t1 - v1.mag2() <= 0) {
241 // std::cerr << "HepLorentzVector::findBoostToCM() - "
242 // << "boostToCM computed for pair of HepLorentzVectors with non-timelike sum"
243 // << std::endl;
244 // // result will make analytic sense but is physically meaningless
245 // }
246  return Hep3Vector(v1 * (-1./t1));
247 } /* boostToCM(w) */
248 
249 } // namespace CLHEP
250 
static const G4double a1
const G4double w[NPOINTSGL]
G4double invariantMass(const G4double E, const ThreeVector &p)
static const double m
Definition: G4SIunits.hh:128