| Geant4 Cross Reference |
1 // -*- C++ -*- 1 // -*- C++ -*-
>> 2 // $Id:$
2 // ------------------------------------------- 3 // ---------------------------------------------------------------------------
3 // 4 //
4 // This file is a part of the CLHEP - a Class 5 // This file is a part of the CLHEP - a Class Library for High Energy Physics.
5 // 6 //
6 // This is the implementation of the HepLorent 7 // This is the implementation of the HepLorentzVector class:
7 // Those methods originating with ZOOM dealing 8 // Those methods originating with ZOOM dealing with comparison (other than
8 // isSpaceLike, isLightlike, isTimelike, which 9 // isSpaceLike, isLightlike, isTimelike, which are in the main part.)
9 // 10 //
10 // 11/29/05 mf in deltaR, replaced the direct 11 // 11/29/05 mf in deltaR, replaced the direct subtraction
11 // pp.phi() - w.getV().phi() with pp.deltaRPhi 12 // pp.phi() - w.getV().phi() with pp.deltaRPhi(w.getV()) which behaves
12 // correctly across the 2pi boundary. 13 // correctly across the 2pi boundary.
13 14
>> 15 #ifdef GNUPRAGMA
>> 16 #pragma implementation
>> 17 #endif
>> 18
14 #include "CLHEP/Vector/LorentzVector.h" 19 #include "CLHEP/Vector/LorentzVector.h"
15 20
16 #include <cmath> 21 #include <cmath>
17 #include <iostream> <<
18 22
19 namespace CLHEP { 23 namespace CLHEP {
20 24
21 //-*********** 25 //-***********
22 // Comparisons 26 // Comparisons
23 //-*********** 27 //-***********
24 28
25 int HepLorentzVector::compare (const HepLorent 29 int HepLorentzVector::compare (const HepLorentzVector & w) const {
26 if ( ee > w.ee ) { 30 if ( ee > w.ee ) {
27 return 1; 31 return 1;
28 } else if ( ee < w.ee ) { 32 } else if ( ee < w.ee ) {
29 return -1; 33 return -1;
30 } else { 34 } else {
31 return ( pp.compare(w.pp) ); 35 return ( pp.compare(w.pp) );
32 } 36 }
33 } /* Compare */ 37 } /* Compare */
34 38
35 bool HepLorentzVector::operator > (const HepLo 39 bool HepLorentzVector::operator > (const HepLorentzVector & w) const {
36 return (compare(w) > 0); 40 return (compare(w) > 0);
37 } 41 }
38 bool HepLorentzVector::operator < (const HepLo 42 bool HepLorentzVector::operator < (const HepLorentzVector & w) const {
39 return (compare(w) < 0); 43 return (compare(w) < 0);
40 } 44 }
41 bool HepLorentzVector::operator>= (const HepLo 45 bool HepLorentzVector::operator>= (const HepLorentzVector & w) const {
42 return (compare(w) >= 0); 46 return (compare(w) >= 0);
43 } 47 }
44 bool HepLorentzVector::operator<= (const HepLo 48 bool HepLorentzVector::operator<= (const HepLorentzVector & w) const {
45 return (compare(w) <= 0); 49 return (compare(w) <= 0);
46 } 50 }
47 51
48 //-******** 52 //-********
49 // isNear 53 // isNear
50 // howNear 54 // howNear
51 //-******** 55 //-********
52 56
53 bool HepLorentzVector::isNear(const HepLorentz 57 bool HepLorentzVector::isNear(const HepLorentzVector & w,
54 double epsilon) const { 58 double epsilon) const {
55 double limit = std::fabs(pp.dot(w.pp)); 59 double limit = std::fabs(pp.dot(w.pp));
56 limit += .25*((ee+w.ee)*(ee+w.ee)); 60 limit += .25*((ee+w.ee)*(ee+w.ee));
57 limit *= epsilon*epsilon; 61 limit *= epsilon*epsilon;
58 double delta = (pp - w.pp).mag2(); 62 double delta = (pp - w.pp).mag2();
59 delta += (ee-w.ee)*(ee-w.ee); 63 delta += (ee-w.ee)*(ee-w.ee);
60 return (delta <= limit ); 64 return (delta <= limit );
61 } /* isNear() */ 65 } /* isNear() */
62 66
63 double HepLorentzVector::howNear(const HepLore 67 double HepLorentzVector::howNear(const HepLorentzVector & w) const {
64 double wdw = std::fabs(pp.dot(w.pp)) + .25*( 68 double wdw = std::fabs(pp.dot(w.pp)) + .25*((ee+w.ee)*(ee+w.ee));
65 double delta = (pp - w.pp).mag2() + (ee-w.ee 69 double delta = (pp - w.pp).mag2() + (ee-w.ee)*(ee-w.ee);
66 if ( (wdw > 0) && (delta < wdw) ) { 70 if ( (wdw > 0) && (delta < wdw) ) {
67 return std::sqrt (delta/wdw); 71 return std::sqrt (delta/wdw);
68 } else if ( (wdw == 0) && (delta == 0) ) { 72 } else if ( (wdw == 0) && (delta == 0) ) {
69 return 0; 73 return 0;
70 } else { 74 } else {
71 return 1; 75 return 1;
72 } 76 }
73 } /* howNear() */ 77 } /* howNear() */
74 78
75 //-********* 79 //-*********
76 // isNearCM 80 // isNearCM
77 // howNearCM 81 // howNearCM
78 //-********* 82 //-*********
79 83
80 bool HepLorentzVector::isNearCM 84 bool HepLorentzVector::isNearCM
81 (const HepLorentzVector & w, double epsi 85 (const HepLorentzVector & w, double epsilon) const {
82 86
83 double tTotal = (ee + w.ee); 87 double tTotal = (ee + w.ee);
84 Hep3Vector vTotal (pp + w.pp); 88 Hep3Vector vTotal (pp + w.pp);
85 double vTotal2 = vTotal.mag2(); 89 double vTotal2 = vTotal.mag2();
86 90
87 if ( vTotal2 >= tTotal*tTotal ) { 91 if ( vTotal2 >= tTotal*tTotal ) {
88 // Either one or both vectors are spacelik 92 // Either one or both vectors are spacelike, or the dominant T components
89 // are in opposite directions. So boostin 93 // are in opposite directions. So boosting and testing makes no sense;
90 // but we do consider two exactly equal ve 94 // but we do consider two exactly equal vectors to be equal in any frame,
91 // even if they are spacelike and can't be 95 // even if they are spacelike and can't be boosted to a CM frame.
92 return (*this == w); 96 return (*this == w);
93 } 97 }
94 98
95 if ( vTotal2 == 0 ) { // no boost needed! 99 if ( vTotal2 == 0 ) { // no boost needed!
96 return (isNear(w, epsilon)); 100 return (isNear(w, epsilon));
97 } 101 }
98 102
99 // Find the boost to the CM frame. We know 103 // Find the boost to the CM frame. We know that the total vector is timelike.
100 104
101 double tRecip = 1./tTotal; 105 double tRecip = 1./tTotal;
102 Hep3Vector bboost ( vTotal * (-tRecip) ); << 106 Hep3Vector boost ( vTotal * (-tRecip) );
103 107
104 //-| Note that you could do pp/t and n 108 //-| Note that you could do pp/t and not be terribly inefficient since
105 //-| SpaceVector/t itself takes 1/t an 109 //-| SpaceVector/t itself takes 1/t and multiplies. The code here saves
106 //-| a redundant check for t=0. 110 //-| a redundant check for t=0.
107 111
108 // Boost both vectors. Since we have the sa 112 // Boost both vectors. Since we have the same boost, there is no need
109 // to repeat the beta and gamma calculation; 113 // to repeat the beta and gamma calculation; and there is no question
110 // about beta >= 1. That is why we don't ju 114 // about beta >= 1. That is why we don't just call w.boosted().
111 115
112 double b2 = vTotal2*tRecip*tRecip; 116 double b2 = vTotal2*tRecip*tRecip;
113 117
114 double ggamma = std::sqrt(1./(1.-b2)); << 118 register double gamma = std::sqrt(1./(1.-b2));
115 double boostDotV1 = bboost.dot(pp); << 119 register double boostDotV1 = boost.dot(pp);
116 double gm1_b2 = (ggamma-1)/b2; << 120 register double gm1_b2 = (gamma-1)/b2;
117 << 121
118 HepLorentzVector w1 ( pp + ((gm1_b2)*boost << 122 HepLorentzVector w1 ( pp + ((gm1_b2)*boostDotV1+gamma*ee) * boost,
119 ggamma * (ee + boostDotV1 << 123 gamma * (ee + boostDotV1) );
120 << 124
121 double boostDotV2 = bboost.dot(w.pp); << 125 register double boostDotV2 = boost.dot(w.pp);
122 HepLorentzVector w2 ( w.pp + ((gm1_b2)*boost << 126 HepLorentzVector w2 ( w.pp + ((gm1_b2)*boostDotV2+gamma*w.ee) * boost,
123 ggamma * (w.ee + boostDot << 127 gamma * (w.ee + boostDotV2) );
124 128
125 return (w1.isNear(w2, epsilon)); 129 return (w1.isNear(w2, epsilon));
126 130
127 } /* isNearCM() */ 131 } /* isNearCM() */
128 132
129 double HepLorentzVector::howNearCM(const HepLo 133 double HepLorentzVector::howNearCM(const HepLorentzVector & w) const {
130 134
131 double tTotal = (ee + w.ee); 135 double tTotal = (ee + w.ee);
132 Hep3Vector vTotal (pp + w.pp); 136 Hep3Vector vTotal (pp + w.pp);
133 double vTotal2 = vTotal.mag2(); 137 double vTotal2 = vTotal.mag2();
134 138
135 if ( vTotal2 >= tTotal*tTotal ) { 139 if ( vTotal2 >= tTotal*tTotal ) {
136 // Either one or both vectors are spacelik 140 // Either one or both vectors are spacelike, or the dominant T components
137 // are in opposite directions. So boostin 141 // are in opposite directions. So boosting and testing makes no sense;
138 // but we do consider two exactly equal ve 142 // but we do consider two exactly equal vectors to be equal in any frame,
139 // even if they are spacelike and can't be 143 // even if they are spacelike and can't be boosted to a CM frame.
140 if (*this == w) { 144 if (*this == w) {
141 return 0; 145 return 0;
142 } else { 146 } else {
143 return 1; 147 return 1;
144 } 148 }
145 } 149 }
146 150
147 if ( vTotal2 == 0 ) { // no boost needed! 151 if ( vTotal2 == 0 ) { // no boost needed!
148 return (howNear(w)); 152 return (howNear(w));
149 } 153 }
150 154
151 // Find the boost to the CM frame. We know 155 // Find the boost to the CM frame. We know that the total vector is timelike.
152 156
153 double tRecip = 1./tTotal; 157 double tRecip = 1./tTotal;
154 Hep3Vector bboost ( vTotal * (-tRecip) ); << 158 Hep3Vector boost ( vTotal * (-tRecip) );
155 159
156 //-| Note that you could do pp/t and n 160 //-| Note that you could do pp/t and not be terribly inefficient since
157 //-| SpaceVector/t itself takes 1/t an 161 //-| SpaceVector/t itself takes 1/t and multiplies. The code here saves
158 //-| a redundant check for t=0. 162 //-| a redundant check for t=0.
159 163
160 // Boost both vectors. Since we have the sa 164 // Boost both vectors. Since we have the same boost, there is no need
161 // to repeat the beta and gamma calculation; 165 // to repeat the beta and gamma calculation; and there is no question
162 // about beta >= 1. That is why we don't ju 166 // about beta >= 1. That is why we don't just call w.boosted().
163 167
164 double b2 = vTotal2*tRecip*tRecip; 168 double b2 = vTotal2*tRecip*tRecip;
165 // if ( b2 >= 1 ) { // NaN-proofing 169 // if ( b2 >= 1 ) { // NaN-proofing
166 // std::cerr << "HepLorentzVector::howNearC 170 // std::cerr << "HepLorentzVector::howNearCM() - "
167 // << "boost vector in howNearCM appears to b 171 // << "boost vector in howNearCM appears to be tachyonic" << std::endl;
168 // } 172 // }
169 double ggamma = std::sqrt(1./(1.-b2)); << 173 register double gamma = std::sqrt(1./(1.-b2));
170 double boostDotV1 = bboost.dot(pp); << 174 register double boostDotV1 = boost.dot(pp);
171 double gm1_b2 = (ggamma-1)/b2; << 175 register double gm1_b2 = (gamma-1)/b2;
172 << 176
173 HepLorentzVector w1 ( pp + ((gm1_b2)*boost << 177 HepLorentzVector w1 ( pp + ((gm1_b2)*boostDotV1+gamma*ee) * boost,
174 ggamma * (ee + boostDotV1 << 178 gamma * (ee + boostDotV1) );
175 << 179
176 double boostDotV2 = bboost.dot(w.pp); << 180 register double boostDotV2 = boost.dot(w.pp);
177 HepLorentzVector w2 ( w.pp + ((gm1_b2)*boost << 181 HepLorentzVector w2 ( w.pp + ((gm1_b2)*boostDotV2+gamma*w.ee) * boost,
178 ggamma * (w.ee + boostDot << 182 gamma * (w.ee + boostDotV2) );
179 183
180 return (w1.howNear(w2)); 184 return (w1.howNear(w2));
181 185
182 } /* howNearCM() */ 186 } /* howNearCM() */
183 187
184 //-************ 188 //-************
185 // deltaR 189 // deltaR
186 // isParallel 190 // isParallel
187 // howParallel 191 // howParallel
188 // howLightlike 192 // howLightlike
189 //-************ 193 //-************
190 194
191 double HepLorentzVector::deltaR ( const HepLor 195 double HepLorentzVector::deltaR ( const HepLorentzVector & w ) const {
192 196
193 double a = eta() - w.eta(); 197 double a = eta() - w.eta();
194 double b = pp.deltaPhi(w.getV()); 198 double b = pp.deltaPhi(w.getV());
195 199
196 return std::sqrt ( a*a + b*b ); 200 return std::sqrt ( a*a + b*b );
197 201
198 } /* deltaR */ 202 } /* deltaR */
199 203
200 // If the difference (in the Euclidean norm) o 204 // If the difference (in the Euclidean norm) of the normalized (in Euclidean
201 // norm) 4-vectors is small, then those 4-vect 205 // norm) 4-vectors is small, then those 4-vectors are considered nearly
202 // parallel. 206 // parallel.
203 207
204 bool HepLorentzVector::isParallel (const HepLo 208 bool HepLorentzVector::isParallel (const HepLorentzVector & w, double epsilon) const {
205 double norm = euclideanNorm(); 209 double norm = euclideanNorm();
206 double wnorm = w.euclideanNorm(); 210 double wnorm = w.euclideanNorm();
207 if ( norm == 0 ) { 211 if ( norm == 0 ) {
208 if ( wnorm == 0 ) { 212 if ( wnorm == 0 ) {
209 return true; 213 return true;
210 } else { 214 } else {
211 return false; 215 return false;
212 } 216 }
213 } 217 }
214 if ( wnorm == 0 ) { 218 if ( wnorm == 0 ) {
215 return false; 219 return false;
216 } 220 }
217 HepLorentzVector w1 = *this / norm; 221 HepLorentzVector w1 = *this / norm;
218 HepLorentzVector w2 = w / wnorm; 222 HepLorentzVector w2 = w / wnorm;
219 return ( (w1-w2).euclideanNorm2() <= epsilon 223 return ( (w1-w2).euclideanNorm2() <= epsilon*epsilon );
220 } /* isParallel */ 224 } /* isParallel */
221 225
222 226
223 double HepLorentzVector::howParallel (const He 227 double HepLorentzVector::howParallel (const HepLorentzVector & w) const {
224 228
225 double norm = euclideanNorm(); 229 double norm = euclideanNorm();
226 double wnorm = w.euclideanNorm(); 230 double wnorm = w.euclideanNorm();
227 if ( norm == 0 ) { 231 if ( norm == 0 ) {
228 if ( wnorm == 0 ) { 232 if ( wnorm == 0 ) {
229 return 0; 233 return 0;
230 } else { 234 } else {
231 return 1; 235 return 1;
232 } 236 }
233 } 237 }
234 if ( wnorm == 0 ) { 238 if ( wnorm == 0 ) {
235 return 1; 239 return 1;
236 } 240 }
237 241
238 HepLorentzVector w1 = *this / norm; 242 HepLorentzVector w1 = *this / norm;
239 HepLorentzVector w2 = w / wnorm; 243 HepLorentzVector w2 = w / wnorm;
240 double x1 = (w1-w2).euclideanNorm(); << 244 double x = (w1-w2).euclideanNorm();
241 return (x1 < 1) ? x1 : 1; << 245 return (x < 1) ? x : 1;
242 246
243 } /* howParallel */ 247 } /* howParallel */
244 248
245 double HepLorentzVector::howLightlike() const 249 double HepLorentzVector::howLightlike() const {
246 double m1 = std::fabs(restMass2()); << 250 double m2 = std::fabs(restMass2());
247 double twoT2 = 2*ee*ee; 251 double twoT2 = 2*ee*ee;
248 if (m1 < twoT2) { << 252 if (m2 < twoT2) {
249 return m1/twoT2; << 253 return m2/twoT2;
250 } else { 254 } else {
251 return 1; 255 return 1;
252 } 256 }
253 } /* HowLightlike */ 257 } /* HowLightlike */
254 258
255 } // namespace CLHEP 259 } // namespace CLHEP
256 260