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