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