CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // -*- C++ -*- 00002 // $Id: LorentzVectorC.cc,v 1.2.4.1 2004/11/30 20:08:40 garren Exp $ 00003 // --------------------------------------------------------------------------- 00004 // 00005 // This file is a part of the CLHEP - a Class Library for High Energy Physics. 00006 // 00007 // This is the implementation of the HepLorentzVector class: 00008 // Those methods originating with ZOOM dealing with comparison (other than 00009 // isSpaceLike, isLightlike, isTimelike, which are in the main part.) 00010 // 00011 // 11/29/05 mf in deltaR, replaced the direct subtraction 00012 // pp.phi() - w.getV().phi() with pp.deltaRPhi(w.getV()) which behaves 00013 // correctly across the 2pi boundary. 00014 00015 #ifdef GNUPRAGMA 00016 #pragma implementation 00017 #endif 00018 00019 #include "CLHEP/Vector/defs.h" 00020 #include "CLHEP/Vector/LorentzVector.h" 00021 00022 #include <cmath> 00023 00024 namespace CLHEP { 00025 00026 //-*********** 00027 // Comparisons 00028 //-*********** 00029 00030 int HepLorentzVector::compare (const HepLorentzVector & w) const { 00031 if ( ee > w.ee ) { 00032 return 1; 00033 } else if ( ee < w.ee ) { 00034 return -1; 00035 } else { 00036 return ( pp.compare(w.pp) ); 00037 } 00038 } /* Compare */ 00039 00040 bool HepLorentzVector::operator > (const HepLorentzVector & w) const { 00041 return (compare(w) > 0); 00042 } 00043 bool HepLorentzVector::operator < (const HepLorentzVector & w) const { 00044 return (compare(w) < 0); 00045 } 00046 bool HepLorentzVector::operator>= (const HepLorentzVector & w) const { 00047 return (compare(w) >= 0); 00048 } 00049 bool HepLorentzVector::operator<= (const HepLorentzVector & w) const { 00050 return (compare(w) <= 0); 00051 } 00052 00053 //-******** 00054 // isNear 00055 // howNear 00056 //-******** 00057 00058 bool HepLorentzVector::isNear(const HepLorentzVector & w, 00059 double epsilon) const { 00060 double limit = fabs(pp.dot(w.pp)); 00061 limit += .25*((ee+w.ee)*(ee+w.ee)); 00062 limit *= epsilon*epsilon; 00063 double delta = (pp - w.pp).mag2(); 00064 delta += (ee-w.ee)*(ee-w.ee); 00065 return (delta <= limit ); 00066 } /* isNear() */ 00067 00068 double HepLorentzVector::howNear(const HepLorentzVector & w) const { 00069 double wdw = fabs(pp.dot(w.pp)) + .25*((ee+w.ee)*(ee+w.ee)); 00070 double delta = (pp - w.pp).mag2() + (ee-w.ee)*(ee-w.ee); 00071 if ( (wdw > 0) && (delta < wdw) ) { 00072 return sqrt (delta/wdw); 00073 } else if ( (wdw == 0) && (delta == 0) ) { 00074 return 0; 00075 } else { 00076 return 1; 00077 } 00078 } /* howNear() */ 00079 00080 //-********* 00081 // isNearCM 00082 // howNearCM 00083 //-********* 00084 00085 bool HepLorentzVector::isNearCM 00086 (const HepLorentzVector & w, double epsilon) const { 00087 00088 double tTotal = (ee + w.ee); 00089 Hep3Vector vTotal (pp + w.pp); 00090 double vTotal2 = vTotal.mag2(); 00091 00092 if ( vTotal2 >= tTotal*tTotal ) { 00093 // Either one or both vectors are spacelike, or the dominant T components 00094 // are in opposite directions. So boosting and testing makes no sense; 00095 // but we do consider two exactly equal vectors to be equal in any frame, 00096 // even if they are spacelike and can't be boosted to a CM frame. 00097 return (*this == w); 00098 } 00099 00100 if ( vTotal2 == 0 ) { // no boost needed! 00101 return (isNear(w, epsilon)); 00102 } 00103 00104 // Find the boost to the CM frame. We know that the total vector is timelike. 00105 00106 double tRecip = 1./tTotal; 00107 Hep3Vector boost ( vTotal * (-tRecip) ); 00108 00109 //-| Note that you could do pp/t and not be terribly inefficient since 00110 //-| SpaceVector/t itself takes 1/t and multiplies. The code here saves 00111 //-| a redundant check for t=0. 00112 00113 // Boost both vectors. Since we have the same boost, there is no need 00114 // to repeat the beta and gamma calculation; and there is no question 00115 // about beta >= 1. That is why we don't just call w.boosted(). 00116 00117 double b2 = vTotal2*tRecip*tRecip; 00118 00119 register double gamma = sqrt(1./(1.-b2)); 00120 register double boostDotV1 = boost.dot(pp); 00121 register double gm1_b2 = (gamma-1)/b2; 00122 00123 HepLorentzVector w1 ( pp + ((gm1_b2)*boostDotV1+gamma*ee) * boost, 00124 gamma * (ee + boostDotV1) ); 00125 00126 register double boostDotV2 = boost.dot(w.pp); 00127 HepLorentzVector w2 ( w.pp + ((gm1_b2)*boostDotV2+gamma*w.ee) * boost, 00128 gamma * (w.ee + boostDotV2) ); 00129 00130 return (w1.isNear(w2, epsilon)); 00131 00132 } /* isNearCM() */ 00133 00134 double HepLorentzVector::howNearCM(const HepLorentzVector & w) const { 00135 00136 double tTotal = (ee + w.ee); 00137 Hep3Vector vTotal (pp + w.pp); 00138 double vTotal2 = vTotal.mag2(); 00139 00140 if ( vTotal2 >= tTotal*tTotal ) { 00141 // Either one or both vectors are spacelike, or the dominant T components 00142 // are in opposite directions. So boosting and testing makes no sense; 00143 // but we do consider two exactly equal vectors to be equal in any frame, 00144 // even if they are spacelike and can't be boosted to a CM frame. 00145 if (*this == w) { 00146 return 0; 00147 } else { 00148 return 1; 00149 } 00150 } 00151 00152 if ( vTotal2 == 0 ) { // no boost needed! 00153 return (howNear(w)); 00154 } 00155 00156 // Find the boost to the CM frame. We know that the total vector is timelike. 00157 00158 double tRecip = 1./tTotal; 00159 Hep3Vector boost ( vTotal * (-tRecip) ); 00160 00161 //-| Note that you could do pp/t and not be terribly inefficient since 00162 //-| SpaceVector/t itself takes 1/t and multiplies. The code here saves 00163 //-| a redundant check for t=0. 00164 00165 // Boost both vectors. Since we have the same boost, there is no need 00166 // to repeat the beta and gamma calculation; and there is no question 00167 // about beta >= 1. That is why we don't just call w.boosted(). 00168 00169 double b2 = vTotal2*tRecip*tRecip; 00170 if ( b2 >= 1 ) { // NaN-proofing 00171 ZMthrowC ( ZMxpvTachyonic ( 00172 "boost vector in howNearCM appears to be tachyonic")); 00173 } 00174 register double gamma = sqrt(1./(1.-b2)); 00175 register double boostDotV1 = boost.dot(pp); 00176 register double gm1_b2 = (gamma-1)/b2; 00177 00178 HepLorentzVector w1 ( pp + ((gm1_b2)*boostDotV1+gamma*ee) * boost, 00179 gamma * (ee + boostDotV1) ); 00180 00181 register double boostDotV2 = boost.dot(w.pp); 00182 HepLorentzVector w2 ( w.pp + ((gm1_b2)*boostDotV2+gamma*w.ee) * boost, 00183 gamma * (w.ee + boostDotV2) ); 00184 00185 return (w1.howNear(w2)); 00186 00187 } /* howNearCM() */ 00188 00189 //-************ 00190 // deltaR 00191 // isParallel 00192 // howParallel 00193 // howLightlike 00194 //-************ 00195 00196 double HepLorentzVector::deltaR ( const HepLorentzVector & w ) const { 00197 00198 double a = eta() - w.eta(); 00199 double b = pp.deltaPhi(w.getV()); 00200 00201 return sqrt ( a*a + b*b ); 00202 00203 } /* deltaR */ 00204 00205 // If the difference (in the Euclidean norm) of the normalized (in Euclidean 00206 // norm) 4-vectors is small, then those 4-vectors are considered nearly 00207 // parallel. 00208 00209 bool HepLorentzVector::isParallel (const HepLorentzVector & w, double epsilon) const { 00210 double norm = euclideanNorm(); 00211 double wnorm = w.euclideanNorm(); 00212 if ( norm == 0 ) { 00213 if ( wnorm == 0 ) { 00214 return true; 00215 } else { 00216 return false; 00217 } 00218 } 00219 if ( wnorm == 0 ) { 00220 return false; 00221 } 00222 HepLorentzVector w1 = *this / norm; 00223 HepLorentzVector w2 = w / wnorm; 00224 return ( (w1-w2).euclideanNorm2() <= epsilon*epsilon ); 00225 } /* isParallel */ 00226 00227 00228 double HepLorentzVector::howParallel (const HepLorentzVector & w) const { 00229 00230 double norm = euclideanNorm(); 00231 double wnorm = w.euclideanNorm(); 00232 if ( norm == 0 ) { 00233 if ( wnorm == 0 ) { 00234 return 0; 00235 } else { 00236 return 1; 00237 } 00238 } 00239 if ( wnorm == 0 ) { 00240 return 1; 00241 } 00242 00243 HepLorentzVector w1 = *this / norm; 00244 HepLorentzVector w2 = w / wnorm; 00245 double x = (w1-w2).euclideanNorm(); 00246 return (x < 1) ? x : 1; 00247 00248 } /* howParallel */ 00249 00250 double HepLorentzVector::howLightlike() const { 00251 double m2 = fabs(restMass2()); 00252 double twoT2 = 2*ee*ee; 00253 if (m2 < twoT2) { 00254 return m2/twoT2; 00255 } else { 00256 return 1; 00257 } 00258 } /* HowLightlike */ 00259 00260 } // namespace CLHEP