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