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