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 
24 namespace 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 
196 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 
209 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 
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
double mag2() const
double deltaR(const HepLorentzVector &v) const
int compare(const HepLorentzVector &w) const
double norm(const HepGenMatrix &m)
Definition: GenMatrix.cc:57
bool operator<(const HepLorentzVector &w) const
int compare(const Hep3Vector &v) const
Definition: SpaceVector.cc:125
double deltaPhi(const Hep3Vector &v2) const
Definition: ThreeVector.cc:172
double dot(const Hep3Vector &) const
#define ZMthrowC(A)
bool operator<=(const HepLorentzVector &w) const
bool isParallel(const HepLorentzVector &w, double epsilon=tolerance) const
bool operator>(const HepLorentzVector &w) const
double euclideanNorm2() const
bool isNearCM(const HepLorentzVector &w, double epsilon=tolerance) const
double mag2() const
double howLightlike() const
double eta() const
Hep3Vector getV() const
double restMass2() const
bool operator>=(const HepLorentzVector &w) const
double howParallel(const HepLorentzVector &w) const
bool isNear(const HepLorentzVector &w, double epsilon=tolerance) const
double howNearCM(const HepLorentzVector &w) const
double howNear(const HepLorentzVector &w) const
double euclideanNorm() const