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