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