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