CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

LorentzVectorC.cc
Go to the documentation of this file.
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
24namespace CLHEP {
25
26//-***********
27// Comparisons
28//-***********
29
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
41 return (compare(w) > 0);
42}
44 return (compare(w) < 0);
45}
47 return (compare(w) >= 0);
48}
50 return (compare(w) <= 0);
51}
52
53//-********
54// isNear
55// howNear
56//-********
57
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
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
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 register double ggamma = std::sqrt(1./(1.-b2));
120 register double boostDotV1 = bboost.dot(pp);
121 register double gm1_b2 = (ggamma-1)/b2;
122
123 HepLorentzVector w1 ( pp + ((gm1_b2)*boostDotV1+ggamma*ee) * bboost,
124 ggamma * (ee + boostDotV1) );
125
126 register 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
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 register double ggamma = std::sqrt(1./(1.-b2));
175 register double boostDotV1 = bboost.dot(pp);
176 register double gm1_b2 = (ggamma-1)/b2;
177
178 HepLorentzVector w1 ( pp + ((gm1_b2)*boostDotV1+ggamma*ee) * bboost,
179 ggamma * (ee + boostDotV1) );
180
181 register 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
196double 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
209bool 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
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
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
#define ZMthrowC(A)
double mag2() const
double dot(const Hep3Vector &) const
double deltaPhi(const Hep3Vector &v2) const
int compare(const Hep3Vector &v) const
bool isParallel(const HepLorentzVector &w, double epsilon=tolerance) const
bool isNearCM(const HepLorentzVector &w, double epsilon=tolerance) const
Hep3Vector getV() const
int compare(const HepLorentzVector &w) const
double howParallel(const HepLorentzVector &w) const
bool operator<=(const HepLorentzVector &w) const
double restMass2() const
double deltaR(const HepLorentzVector &v) const
double euclideanNorm() const
double howNear(const HepLorentzVector &w) const
double howNearCM(const HepLorentzVector &w) const
bool operator<(const HepLorentzVector &w) const
double euclideanNorm2() const
bool isNear(const HepLorentzVector &w, double epsilon=tolerance) const
bool operator>(const HepLorentzVector &w) const
bool operator>=(const HepLorentzVector &w) const
double norm(const HepGenMatrix &m)
Definition GenMatrix.cc:57