Geant4 Cross Reference |
1 // -*- C++ -*- 1 2 // ------------------------------------------- 3 // 4 // This file is a part of the CLHEP - a Class 5 // 6 // This is the implementation of the HepLorent 7 // Those methods originating with ZOOM dealing 8 // isSpaceLike, isLightlike, isTimelike, which 9 // 10 // 11/29/05 mf in deltaR, replaced the direct 11 // pp.phi() - w.getV().phi() with pp.deltaRPhi 12 // correctly across the 2pi boundary. 13 14 #include "CLHEP/Vector/LorentzVector.h" 15 16 #include <cmath> 17 #include <iostream> 18 19 namespace CLHEP { 20 21 //-*********** 22 // Comparisons 23 //-*********** 24 25 int HepLorentzVector::compare (const HepLorent 26 if ( ee > w.ee ) { 27 return 1; 28 } else if ( ee < w.ee ) { 29 return -1; 30 } else { 31 return ( pp.compare(w.pp) ); 32 } 33 } /* Compare */ 34 35 bool HepLorentzVector::operator > (const HepLo 36 return (compare(w) > 0); 37 } 38 bool HepLorentzVector::operator < (const HepLo 39 return (compare(w) < 0); 40 } 41 bool HepLorentzVector::operator>= (const HepLo 42 return (compare(w) >= 0); 43 } 44 bool HepLorentzVector::operator<= (const HepLo 45 return (compare(w) <= 0); 46 } 47 48 //-******** 49 // isNear 50 // howNear 51 //-******** 52 53 bool HepLorentzVector::isNear(const HepLorentz 54 double epsilon) const { 55 double limit = std::fabs(pp.dot(w.pp)); 56 limit += .25*((ee+w.ee)*(ee+w.ee)); 57 limit *= epsilon*epsilon; 58 double delta = (pp - w.pp).mag2(); 59 delta += (ee-w.ee)*(ee-w.ee); 60 return (delta <= limit ); 61 } /* isNear() */ 62 63 double HepLorentzVector::howNear(const HepLore 64 double wdw = std::fabs(pp.dot(w.pp)) + .25*( 65 double delta = (pp - w.pp).mag2() + (ee-w.ee 66 if ( (wdw > 0) && (delta < wdw) ) { 67 return std::sqrt (delta/wdw); 68 } else if ( (wdw == 0) && (delta == 0) ) { 69 return 0; 70 } else { 71 return 1; 72 } 73 } /* howNear() */ 74 75 //-********* 76 // isNearCM 77 // howNearCM 78 //-********* 79 80 bool HepLorentzVector::isNearCM 81 (const HepLorentzVector & w, double epsi 82 83 double tTotal = (ee + w.ee); 84 Hep3Vector vTotal (pp + w.pp); 85 double vTotal2 = vTotal.mag2(); 86 87 if ( vTotal2 >= tTotal*tTotal ) { 88 // Either one or both vectors are spacelik 89 // are in opposite directions. So boostin 90 // but we do consider two exactly equal ve 91 // even if they are spacelike and can't be 92 return (*this == w); 93 } 94 95 if ( vTotal2 == 0 ) { // no boost needed! 96 return (isNear(w, epsilon)); 97 } 98 99 // Find the boost to the CM frame. We know 100 101 double tRecip = 1./tTotal; 102 Hep3Vector bboost ( vTotal * (-tRecip) ); 103 104 //-| Note that you could do pp/t and n 105 //-| SpaceVector/t itself takes 1/t an 106 //-| a redundant check for t=0. 107 108 // Boost both vectors. Since we have the sa 109 // to repeat the beta and gamma calculation; 110 // about beta >= 1. That is why we don't ju 111 112 double b2 = vTotal2*tRecip*tRecip; 113 114 double ggamma = std::sqrt(1./(1.-b2)); 115 double boostDotV1 = bboost.dot(pp); 116 double gm1_b2 = (ggamma-1)/b2; 117 118 HepLorentzVector w1 ( pp + ((gm1_b2)*boost 119 ggamma * (ee + boostDotV1 120 121 double boostDotV2 = bboost.dot(w.pp); 122 HepLorentzVector w2 ( w.pp + ((gm1_b2)*boost 123 ggamma * (w.ee + boostDot 124 125 return (w1.isNear(w2, epsilon)); 126 127 } /* isNearCM() */ 128 129 double HepLorentzVector::howNearCM(const HepLo 130 131 double tTotal = (ee + w.ee); 132 Hep3Vector vTotal (pp + w.pp); 133 double vTotal2 = vTotal.mag2(); 134 135 if ( vTotal2 >= tTotal*tTotal ) { 136 // Either one or both vectors are spacelik 137 // are in opposite directions. So boostin 138 // but we do consider two exactly equal ve 139 // even if they are spacelike and can't be 140 if (*this == w) { 141 return 0; 142 } else { 143 return 1; 144 } 145 } 146 147 if ( vTotal2 == 0 ) { // no boost needed! 148 return (howNear(w)); 149 } 150 151 // Find the boost to the CM frame. We know 152 153 double tRecip = 1./tTotal; 154 Hep3Vector bboost ( vTotal * (-tRecip) ); 155 156 //-| Note that you could do pp/t and n 157 //-| SpaceVector/t itself takes 1/t an 158 //-| a redundant check for t=0. 159 160 // Boost both vectors. Since we have the sa 161 // to repeat the beta and gamma calculation; 162 // about beta >= 1. That is why we don't ju 163 164 double b2 = vTotal2*tRecip*tRecip; 165 // if ( b2 >= 1 ) { // NaN-proofing 166 // std::cerr << "HepLorentzVector::howNearC 167 // << "boost vector in howNearCM appears to b 168 // } 169 double ggamma = std::sqrt(1./(1.-b2)); 170 double boostDotV1 = bboost.dot(pp); 171 double gm1_b2 = (ggamma-1)/b2; 172 173 HepLorentzVector w1 ( pp + ((gm1_b2)*boost 174 ggamma * (ee + boostDotV1 175 176 double boostDotV2 = bboost.dot(w.pp); 177 HepLorentzVector w2 ( w.pp + ((gm1_b2)*boost 178 ggamma * (w.ee + boostDot 179 180 return (w1.howNear(w2)); 181 182 } /* howNearCM() */ 183 184 //-************ 185 // deltaR 186 // isParallel 187 // howParallel 188 // howLightlike 189 //-************ 190 191 double HepLorentzVector::deltaR ( const HepLor 192 193 double a = eta() - w.eta(); 194 double b = pp.deltaPhi(w.getV()); 195 196 return std::sqrt ( a*a + b*b ); 197 198 } /* deltaR */ 199 200 // If the difference (in the Euclidean norm) o 201 // norm) 4-vectors is small, then those 4-vect 202 // parallel. 203 204 bool HepLorentzVector::isParallel (const HepLo 205 double norm = euclideanNorm(); 206 double wnorm = w.euclideanNorm(); 207 if ( norm == 0 ) { 208 if ( wnorm == 0 ) { 209 return true; 210 } else { 211 return false; 212 } 213 } 214 if ( wnorm == 0 ) { 215 return false; 216 } 217 HepLorentzVector w1 = *this / norm; 218 HepLorentzVector w2 = w / wnorm; 219 return ( (w1-w2).euclideanNorm2() <= epsilon 220 } /* isParallel */ 221 222 223 double HepLorentzVector::howParallel (const He 224 225 double norm = euclideanNorm(); 226 double wnorm = w.euclideanNorm(); 227 if ( norm == 0 ) { 228 if ( wnorm == 0 ) { 229 return 0; 230 } else { 231 return 1; 232 } 233 } 234 if ( wnorm == 0 ) { 235 return 1; 236 } 237 238 HepLorentzVector w1 = *this / norm; 239 HepLorentzVector w2 = w / wnorm; 240 double x1 = (w1-w2).euclideanNorm(); 241 return (x1 < 1) ? x1 : 1; 242 243 } /* howParallel */ 244 245 double HepLorentzVector::howLightlike() const 246 double m1 = std::fabs(restMass2()); 247 double twoT2 = 2*ee*ee; 248 if (m1 < twoT2) { 249 return m1/twoT2; 250 } else { 251 return 1; 252 } 253 } /* HowLightlike */ 254 255 } // namespace CLHEP 256