Geant4 Cross Reference |
1 // -*- C++ -*- 1 // -*- C++ -*- 2 // ------------------------------------------- 2 // --------------------------------------------------------------------------- 3 // 3 // 4 // This file is a part of the CLHEP - a Class 4 // This file is a part of the CLHEP - a Class Library for High Energy Physics. 5 // 5 // 6 // This is part of the implementation of the H 6 // This is part of the implementation of the HepLorentzVector class: 7 // Those methods which originated from ZOOM an 7 // Those methods which originated from ZOOM and which deal with relativistic 8 // kinematic properties. 8 // kinematic properties. 9 // 9 // 10 10 >> 11 #ifdef GNUPRAGMA >> 12 #pragma implementation >> 13 #endif >> 14 11 #include "CLHEP/Vector/LorentzVector.h" 15 #include "CLHEP/Vector/LorentzVector.h" 12 16 13 #include <cmath> 17 #include <cmath> 14 #include <iostream> << 15 18 16 namespace CLHEP { 19 namespace CLHEP { 17 20 18 //-****************** 21 //-****************** 19 // Metric flexibility 22 // Metric flexibility 20 //-****************** 23 //-****************** 21 24 22 ZMpvMetric_t HepLorentzVector::setMetric( ZMpv 25 ZMpvMetric_t HepLorentzVector::setMetric( ZMpvMetric_t a1 ) { 23 ZMpvMetric_t oldMetric = (metric > 0) ? Time 26 ZMpvMetric_t oldMetric = (metric > 0) ? TimePositive : TimeNegative; 24 if ( a1 == TimeNegative ) { 27 if ( a1 == TimeNegative ) { 25 metric = -1.0; 28 metric = -1.0; 26 } else { 29 } else { 27 metric = 1.0; 30 metric = 1.0; 28 } 31 } 29 return oldMetric; 32 return oldMetric; 30 } 33 } 31 34 32 ZMpvMetric_t HepLorentzVector::getMetric() { 35 ZMpvMetric_t HepLorentzVector::getMetric() { 33 return ( (metric > 0) ? TimePositive : TimeN 36 return ( (metric > 0) ? TimePositive : TimeNegative ); 34 } 37 } 35 38 36 //-******** 39 //-******** 37 // plus 40 // plus 38 // minus 41 // minus 39 //-******** 42 //-******** 40 43 41 double HepLorentzVector::plus (const Hep3Vecto 44 double HepLorentzVector::plus (const Hep3Vector & ref) const { 42 double r = ref.mag(); 45 double r = ref.mag(); 43 if (r == 0) { 46 if (r == 0) { 44 std::cerr << "HepLorentzVector::plus() - " 47 std::cerr << "HepLorentzVector::plus() - " 45 << "A zero vector used as reference to L 48 << "A zero vector used as reference to LorentzVector plus-part" 46 << std::endl; 49 << std::endl; 47 return ee; 50 return ee; 48 } 51 } 49 return ee + pp.dot(ref)/r; 52 return ee + pp.dot(ref)/r; 50 } /* plus */ 53 } /* plus */ 51 54 52 double HepLorentzVector::minus (const Hep3Vect 55 double HepLorentzVector::minus (const Hep3Vector & ref) const { 53 double r = ref.mag(); 56 double r = ref.mag(); 54 if (r == 0) { 57 if (r == 0) { 55 std::cerr << "HepLorentzVector::minus() - 58 std::cerr << "HepLorentzVector::minus() - " 56 << "A zero vector used as reference to L 59 << "A zero vector used as reference to LorentzVector minus-part" 57 << std::endl; 60 << std::endl; 58 return ee; 61 return ee; 59 } 62 } 60 return ee - pp.dot(ref)/r; 63 return ee - pp.dot(ref)/r; 61 } /* plus */ 64 } /* plus */ 62 65 63 HepLorentzVector HepLorentzVector::rest4Vector 66 HepLorentzVector HepLorentzVector::rest4Vector() const { 64 return HepLorentzVector (0, 0, 0, (t() < 0.0 67 return HepLorentzVector (0, 0, 0, (t() < 0.0 ? -m() : m())); 65 } 68 } 66 69 67 //-******** 70 //-******** 68 // beta 71 // beta 69 // gamma 72 // gamma 70 //-******** 73 //-******** 71 74 72 double HepLorentzVector::beta() const { 75 double HepLorentzVector::beta() const { 73 if (ee == 0) { 76 if (ee == 0) { 74 if (pp.mag2() == 0) { 77 if (pp.mag2() == 0) { 75 return 0; 78 return 0; 76 } else { 79 } else { 77 std::cerr << "HepLorentzVector::beta() - 80 std::cerr << "HepLorentzVector::beta() - " 78 << "beta computed for HepLorentzVector 81 << "beta computed for HepLorentzVector with t=0 -- infinite result" 79 << std::endl; 82 << std::endl; 80 return 1./ee; 83 return 1./ee; 81 } 84 } 82 } 85 } 83 // if (restMass2() <= 0) { 86 // if (restMass2() <= 0) { 84 // std::cerr << "HepLorentzVector::beta() - 87 // std::cerr << "HepLorentzVector::beta() - " 85 // << "beta computed for a non-timelike H 88 // << "beta computed for a non-timelike HepLorentzVector" << std::endl; 86 // // result will make analytic sense b 89 // // result will make analytic sense but is physically meaningless 87 // } 90 // } 88 return std::sqrt (pp.mag2() / (ee*ee)) ; 91 return std::sqrt (pp.mag2() / (ee*ee)) ; 89 } /* beta */ 92 } /* beta */ 90 93 91 double HepLorentzVector::gamma() const { 94 double HepLorentzVector::gamma() const { 92 double v2 = pp.mag2(); 95 double v2 = pp.mag2(); 93 double t2 = ee*ee; 96 double t2 = ee*ee; 94 if (ee == 0) { 97 if (ee == 0) { 95 if (pp.mag2() == 0) { 98 if (pp.mag2() == 0) { 96 return 1; 99 return 1; 97 } else { 100 } else { 98 std::cerr << "HepLorentzVector::gamma() 101 std::cerr << "HepLorentzVector::gamma() - " 99 << "gamma computed for HepLorentzVecto 102 << "gamma computed for HepLorentzVector with t=0 -- zero result" 100 << std::endl; 103 << std::endl; 101 return 0; 104 return 0; 102 } 105 } 103 } 106 } 104 if (t2 < v2) { 107 if (t2 < v2) { 105 std::cerr << "HepLorentzVector::gamma() - 108 std::cerr << "HepLorentzVector::gamma() - " 106 << "gamma computed for a spacelike HepLo 109 << "gamma computed for a spacelike HepLorentzVector -- imaginary result" 107 << std::endl; 110 << std::endl; 108 // analytic result would be imaginary. 111 // analytic result would be imaginary. 109 return 0; 112 return 0; 110 // } else if ( t2 == v2 ) { 113 // } else if ( t2 == v2 ) { 111 // std::cerr << "HepLorentzVector::gamma() 114 // std::cerr << "HepLorentzVector::gamma() - " 112 // << "gamma computed for a lightlike Hep 115 // << "gamma computed for a lightlike HepLorentzVector -- infinite result" 113 // << std::endl; 116 // << std::endl; 114 } 117 } 115 return 1./std::sqrt(1. - v2/t2 ); 118 return 1./std::sqrt(1. - v2/t2 ); 116 } /* gamma */ 119 } /* gamma */ 117 120 118 121 119 //-*************** 122 //-*************** 120 // rapidity 123 // rapidity 121 // pseudorapidity 124 // pseudorapidity 122 // eta 125 // eta 123 //-*************** 126 //-*************** 124 127 125 double HepLorentzVector::rapidity() const { 128 double HepLorentzVector::rapidity() const { 126 double z1 = pp.getZ(); 129 double z1 = pp.getZ(); 127 // if (std::fabs(ee) == std::fabs(z1)) { 130 // if (std::fabs(ee) == std::fabs(z1)) { 128 // std::cerr << "HepLorentzVector::rapidity 131 // std::cerr << "HepLorentzVector::rapidity() - " 129 // << "rapidity for 4-vector with |E| = | 132 // << "rapidity for 4-vector with |E| = |Pz| -- infinite result" 130 // << std::endl; 133 // << std::endl; 131 // } 134 // } 132 if (std::fabs(ee) < std::fabs(z1)) { 135 if (std::fabs(ee) < std::fabs(z1)) { 133 std::cerr << "HepLorentzVector::rapidity() 136 std::cerr << "HepLorentzVector::rapidity() - " 134 << "rapidity for spacelike 4-vector with 137 << "rapidity for spacelike 4-vector with |E| < |Pz| -- undefined" 135 << std::endl; 138 << std::endl; 136 return 0; 139 return 0; 137 } 140 } 138 double q = (ee + z1) / (ee - z1); 141 double q = (ee + z1) / (ee - z1); 139 //-| This cannot be negative now, sinc 142 //-| This cannot be negative now, since both numerator 140 //-| and denominator have the same sig 143 //-| and denominator have the same sign as ee. 141 return .5 * std::log(q); 144 return .5 * std::log(q); 142 } /* rapidity */ 145 } /* rapidity */ 143 146 144 double HepLorentzVector::rapidity(const Hep3Ve 147 double HepLorentzVector::rapidity(const Hep3Vector & ref) const { 145 double r = ref.mag2(); 148 double r = ref.mag2(); 146 if (r == 0) { 149 if (r == 0) { 147 std::cerr << "HepLorentzVector::rapidity() 150 std::cerr << "HepLorentzVector::rapidity() - " 148 << "A zero vector used as reference to L 151 << "A zero vector used as reference to LorentzVector rapidity" 149 << std::endl; 152 << std::endl; 150 return 0; 153 return 0; 151 } 154 } 152 double vdotu = pp.dot(ref)/std::sqrt(r); 155 double vdotu = pp.dot(ref)/std::sqrt(r); 153 // if (std::fabs(ee) == std::fabs(vdotu)) { 156 // if (std::fabs(ee) == std::fabs(vdotu)) { 154 // std::cerr << "HepLorentzVector::rapidity 157 // std::cerr << "HepLorentzVector::rapidity() - " 155 // << "rapidity for 4-vector with |E| = | 158 // << "rapidity for 4-vector with |E| = |Pu| -- infinite result" 156 // << std::endl; 159 // << std::endl; 157 // } 160 // } 158 if (std::fabs(ee) < std::fabs(vdotu)) { 161 if (std::fabs(ee) < std::fabs(vdotu)) { 159 std::cerr << "HepLorentzVector::rapidity() 162 std::cerr << "HepLorentzVector::rapidity() - " 160 << "rapidity for spacelike 4-vector with 163 << "rapidity for spacelike 4-vector with |E| < |P*ref| -- undefined " 161 << std::endl; 164 << std::endl; 162 return 0; 165 return 0; 163 } 166 } 164 double q = (ee + vdotu) / (ee - vdotu); 167 double q = (ee + vdotu) / (ee - vdotu); 165 return .5 * std::log(q); 168 return .5 * std::log(q); 166 } /* rapidity(ref) */ 169 } /* rapidity(ref) */ 167 170 168 double HepLorentzVector::coLinearRapidity() co 171 double HepLorentzVector::coLinearRapidity() const { 169 double v1 = pp.mag(); 172 double v1 = pp.mag(); 170 // if (std::fabs(ee) == std::fabs(v1)) { 173 // if (std::fabs(ee) == std::fabs(v1)) { 171 // std::cerr << "HepLorentzVector::coLinear 174 // std::cerr << "HepLorentzVector::coLinearRapidity() - " 172 // << "co-Linear rapidity for 4-vector wi 175 // << "co-Linear rapidity for 4-vector with |E| = |P| -- infinite result" 173 // << std::endl; 176 // << std::endl; 174 // } 177 // } 175 if (std::fabs(ee) < std::fabs(v1)) { 178 if (std::fabs(ee) < std::fabs(v1)) { 176 std::cerr << "HepLorentzVector::coLinearRa 179 std::cerr << "HepLorentzVector::coLinearRapidity() - " 177 << "co-linear rapidity for spacelike 4-v 180 << "co-linear rapidity for spacelike 4-vector -- undefined" 178 << std::endl; 181 << std::endl; 179 return 0; 182 return 0; 180 } 183 } 181 double q = (ee + v1) / (ee - v1); 184 double q = (ee + v1) / (ee - v1); 182 return .5 * std::log(q); 185 return .5 * std::log(q); 183 } /* rapidity */ 186 } /* rapidity */ 184 187 185 //-************* 188 //-************* 186 // invariantMass 189 // invariantMass 187 //-************* 190 //-************* 188 191 189 double HepLorentzVector::invariantMass(const H 192 double HepLorentzVector::invariantMass(const HepLorentzVector & w) const { 190 double m1 = invariantMass2(w); 193 double m1 = invariantMass2(w); 191 if (m1 < 0) { 194 if (m1 < 0) { 192 // We should find out why: 195 // We should find out why: 193 if ( ee * w.ee < 0 ) { 196 if ( ee * w.ee < 0 ) { 194 std::cerr << "HepLorentzVector::invarian 197 std::cerr << "HepLorentzVector::invariantMass() - " 195 << "invariant mass meaningless: \n" 198 << "invariant mass meaningless: \n" 196 << "a negative-mass input led to space 199 << "a negative-mass input led to spacelike 4-vector sum" << std::endl; 197 return 0; 200 return 0; 198 } else if ( (isSpacelike() && !isLightlike 201 } else if ( (isSpacelike() && !isLightlike()) || 199 (w.isSpacelike() && !w.isLight 202 (w.isSpacelike() && !w.isLightlike()) ) { 200 std::cerr << "HepLorentzVector::invari 203 std::cerr << "HepLorentzVector::invariantMass() - " 201 << "invariant mass meaningless becau 204 << "invariant mass meaningless because of spacelike input" 202 << std::endl; 205 << std::endl; 203 return 0; 206 return 0; 204 } else { 207 } else { 205 // Invariant mass squared for a pair of 208 // Invariant mass squared for a pair of timelike or lightlike vectors 206 // mathematically cannot be negative. I 209 // mathematically cannot be negative. If the vectors are within the 207 // tolerance of being lightlike or timel 210 // tolerance of being lightlike or timelike, we can assume that prior 208 // or current roundoffs have caused the 211 // or current roundoffs have caused the negative result, and return 0 209 // without comment. 212 // without comment. 210 return 0; 213 return 0; 211 } 214 } 212 } 215 } 213 return (ee+w.ee >=0 ) ? std::sqrt(m1) : - st 216 return (ee+w.ee >=0 ) ? std::sqrt(m1) : - std::sqrt(m1); 214 } /* invariantMass */ 217 } /* invariantMass */ 215 218 216 //-*************** 219 //-*************** 217 // findBoostToCM 220 // findBoostToCM 218 //-*************** 221 //-*************** 219 222 220 Hep3Vector HepLorentzVector::findBoostToCM() c 223 Hep3Vector HepLorentzVector::findBoostToCM() const { 221 return -boostVector(); 224 return -boostVector(); 222 } /* boostToCM() */ 225 } /* boostToCM() */ 223 226 224 Hep3Vector HepLorentzVector::findBoostToCM (co 227 Hep3Vector HepLorentzVector::findBoostToCM (const HepLorentzVector & w) const { 225 double t1 = ee + w.ee; 228 double t1 = ee + w.ee; 226 Hep3Vector v1 = pp + w.pp; 229 Hep3Vector v1 = pp + w.pp; 227 if (t1 == 0) { 230 if (t1 == 0) { 228 if (v1.mag2() == 0) { 231 if (v1.mag2() == 0) { 229 return Hep3Vector(0,0,0); 232 return Hep3Vector(0,0,0); 230 } else { 233 } else { 231 std::cerr << "HepLorentzVector::findBoos 234 std::cerr << "HepLorentzVector::findBoostToCM() - " 232 << "boostToCM computed for two 4-vecto 235 << "boostToCM computed for two 4-vectors with combined t=0 -- " 233 << "infinite result" << std::endl; 236 << "infinite result" << std::endl; 234 return Hep3Vector(v1*(1./t1)); // Yup, 1 237 return Hep3Vector(v1*(1./t1)); // Yup, 1/0 -- that is how we return infinity 235 } 238 } 236 } 239 } 237 // if (t1*t1 - v1.mag2() <= 0) { 240 // if (t1*t1 - v1.mag2() <= 0) { 238 // std::cerr << "HepLorentzVector::findBoos 241 // std::cerr << "HepLorentzVector::findBoostToCM() - " 239 // << "boostToCM computed for pair of He 242 // << "boostToCM computed for pair of HepLorentzVectors with non-timelike sum" 240 // << std::endl; 243 // << std::endl; 241 // // result will make analytic sense b 244 // // result will make analytic sense but is physically meaningless 242 // } 245 // } 243 return Hep3Vector(v1 * (-1./t1)); 246 return Hep3Vector(v1 * (-1./t1)); 244 } /* boostToCM(w) */ 247 } /* boostToCM(w) */ 245 248 246 } // namespace CLHEP 249 } // namespace CLHEP 247 250 248 251