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