00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifdef GNUPRAGMA
00013 #pragma implementation
00014 #endif
00015
00016 #include "CLHEP/Vector/LorentzVector.h"
00017
00018 #include <iostream>
00019
00020 namespace CLHEP {
00021
00022 double HepLorentzVector::tolerance =
00023 Hep3Vector::ToleranceTicks * 2.22045e-16;
00024 double HepLorentzVector::metric = 1.0;
00025
00026 double HepLorentzVector::operator () (int i) const {
00027 switch(i) {
00028 case X:
00029 case Y:
00030 case Z:
00031 return pp(i);
00032 case T:
00033 return e();
00034 default:
00035 std::cerr << "HepLorentzVector subscripting: bad index (" << i << ")"
00036 << std::endl;
00037 }
00038 return 0.;
00039 }
00040
00041 double & HepLorentzVector::operator () (int i) {
00042 static double dummy;
00043 switch(i) {
00044 case X:
00045 case Y:
00046 case Z:
00047 return pp(i);
00048 case T:
00049 return ee;
00050 default:
00051 std::cerr
00052 << "HepLorentzVector subscripting: bad index (" << i << ")"
00053 << std::endl;
00054 return dummy;
00055 }
00056 }
00057
00058 HepLorentzVector & HepLorentzVector::boost
00059 (double bx, double by, double bz){
00060 double b2 = bx*bx + by*by + bz*bz;
00061 register double ggamma = 1.0 / std::sqrt(1.0 - b2);
00062 register double bp = bx*x() + by*y() + bz*z();
00063 register double gamma2 = b2 > 0 ? (ggamma - 1.0)/b2 : 0.0;
00064
00065 setX(x() + gamma2*bp*bx + ggamma*bx*t());
00066 setY(y() + gamma2*bp*by + ggamma*by*t());
00067 setZ(z() + gamma2*bp*bz + ggamma*bz*t());
00068 setT(ggamma*(t() + bp));
00069 return *this;
00070 }
00071
00072 HepLorentzVector & HepLorentzVector::rotateX(double a) {
00073 pp.rotateX(a);
00074 return *this;
00075 }
00076 HepLorentzVector & HepLorentzVector::rotateY(double a) {
00077 pp.rotateY(a);
00078 return *this;
00079 }
00080 HepLorentzVector & HepLorentzVector::rotateZ(double a) {
00081 pp.rotateZ(a);
00082 return *this;
00083 }
00084
00085 HepLorentzVector & HepLorentzVector::rotateUz(const Hep3Vector &v1) {
00086 pp.rotateUz(v1);
00087 return *this;
00088 }
00089
00090 std::ostream & operator<< (std::ostream & os, const HepLorentzVector & v1)
00091 {
00092 return os << "(" << v1.x() << "," << v1.y() << "," << v1.z()
00093 << ";" << v1.t() << ")";
00094 }
00095
00096 std::istream & operator>> (std::istream & is, HepLorentzVector & v1) {
00097
00098
00099
00100
00101
00102
00103 double x, y, z, t;
00104 char c;
00105
00106 is >> std::ws >> c;
00107
00108
00109 if (is.fail() || c != '(' ) {
00110 std::cerr << "Could not find required opening parenthesis "
00111 << "in input of a HepLorentzVector" << std::endl;
00112 return is;
00113 }
00114
00115 is >> x >> std::ws >> c;
00116 if (is.fail() || c != ',' ) {
00117 std::cerr << "Could not find x value and required trailing comma "
00118 << "in input of a HepLorentzVector" << std::endl;
00119 return is;
00120 }
00121
00122 is >> y >> std::ws >> c;
00123 if (is.fail() || c != ',' ) {
00124 std::cerr << "Could not find y value and required trailing comma "
00125 << "in input of a HepLorentzVector" << std::endl;
00126 return is;
00127 }
00128
00129 is >> z >> std::ws >> c;
00130 if (is.fail() || c != ';' ) {
00131 std::cerr << "Could not find z value and required trailing semicolon "
00132 << "in input of a HepLorentzVector" << std::endl;
00133 return is;
00134 }
00135
00136 is >> t >> std::ws >> c;
00137 if (is.fail() || c != ')' ) {
00138 std::cerr << "Could not find t value and required close parenthesis "
00139 << "in input of a HepLorentzVector" << std::endl;
00140 return is;
00141 }
00142
00143 v1.setX(x);
00144 v1.setY(y);
00145 v1.setZ(z);
00146 v1.setT(t);
00147 return is;
00148 }
00149
00150
00151
00152 HepLorentzVector & HepLorentzVector::operator /= (double c) {
00153
00154
00155
00156
00157
00158
00159 double oneOverC = 1.0/c;
00160 pp *= oneOverC;
00161 ee *= oneOverC;
00162 return *this;
00163 }
00164
00165 HepLorentzVector operator / (const HepLorentzVector & w, double c) {
00166
00167
00168
00169
00170
00171
00172 double oneOverC = 1.0/c;
00173 return HepLorentzVector (w.getV() * oneOverC,
00174 w.getT() * oneOverC);
00175 }
00176
00177 Hep3Vector HepLorentzVector::boostVector() const {
00178 if (ee == 0) {
00179 if (pp.mag2() == 0) {
00180 return Hep3Vector(0,0,0);
00181 } else {
00182 std::cerr << "HepLorentzVector::boostVector() - "
00183 << "boostVector computed for LorentzVector with t=0 -- infinite result"
00184 << std::endl;
00185 return pp/ee;
00186 }
00187 }
00188 if (restMass2() <= 0) {
00189 std::cerr << "HepLorentzVector::boostVector() - "
00190 << "boostVector computed for a non-timelike LorentzVector " << std::endl;
00191
00192 }
00193 return pp * (1./ee);
00194 }
00195
00196
00197 HepLorentzVector & HepLorentzVector::boostX (double bbeta){
00198 register double b2 = bbeta*bbeta;
00199 if (b2 >= 1) {
00200 std::cerr << "HepLorentzVector::boostX() - "
00201 << "boost along X with beta >= 1 (speed of light) -- \n"
00202 << "no boost done" << std::endl;
00203 } else {
00204 register double ggamma = std::sqrt(1./(1-b2));
00205 register double tt = ee;
00206 ee = ggamma*(ee + bbeta*pp.getX());
00207 pp.setX(ggamma*(pp.getX() + bbeta*tt));
00208 }
00209 return *this;
00210 }
00211
00212 HepLorentzVector & HepLorentzVector::boostY (double bbeta){
00213 register double b2 = bbeta*bbeta;
00214 if (b2 >= 1) {
00215 std::cerr << "HepLorentzVector::boostY() - "
00216 << "boost along Y with beta >= 1 (speed of light) -- \n"
00217 << "no boost done" << std::endl;
00218 } else {
00219 register double ggamma = std::sqrt(1./(1-b2));
00220 register double tt = ee;
00221 ee = ggamma*(ee + bbeta*pp.getY());
00222 pp.setY(ggamma*(pp.getY() + bbeta*tt));
00223 }
00224 return *this;
00225 }
00226
00227 HepLorentzVector & HepLorentzVector::boostZ (double bbeta){
00228 register double b2 = bbeta*bbeta;
00229 if (b2 >= 1) {
00230 std::cerr << "HepLorentzVector::boostZ() - "
00231 << "boost along Z with beta >= 1 (speed of light) -- \n"
00232 << "no boost done" << std::endl;
00233 } else {
00234 register double ggamma = std::sqrt(1./(1-b2));
00235 register double tt = ee;
00236 ee = ggamma*(ee + bbeta*pp.getZ());
00237 pp.setZ(ggamma*(pp.getZ() + bbeta*tt));
00238 }
00239 return *this;
00240 }
00241
00242 double HepLorentzVector::setTolerance ( double tol ) {
00243
00244 double oldTolerance (tolerance);
00245 tolerance = tol;
00246 return oldTolerance;
00247 }
00248
00249 double HepLorentzVector::getTolerance ( ) {
00250
00251 return tolerance;
00252 }
00253
00254 }