00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015 #ifdef GNUPRAGMA
00016 #pragma implementation
00017 #endif
00018
00019 #include "CLHEP/Vector/LorentzVector.h"
00020
00021 #include <cmath>
00022
00023 namespace CLHEP {
00024
00025
00026
00027
00028
00029 int HepLorentzVector::compare (const HepLorentzVector & w) const {
00030 if ( ee > w.ee ) {
00031 return 1;
00032 } else if ( ee < w.ee ) {
00033 return -1;
00034 } else {
00035 return ( pp.compare(w.pp) );
00036 }
00037 }
00038
00039 bool HepLorentzVector::operator > (const HepLorentzVector & w) const {
00040 return (compare(w) > 0);
00041 }
00042 bool HepLorentzVector::operator < (const HepLorentzVector & w) const {
00043 return (compare(w) < 0);
00044 }
00045 bool HepLorentzVector::operator>= (const HepLorentzVector & w) const {
00046 return (compare(w) >= 0);
00047 }
00048 bool HepLorentzVector::operator<= (const HepLorentzVector & w) const {
00049 return (compare(w) <= 0);
00050 }
00051
00052
00053
00054
00055
00056
00057 bool HepLorentzVector::isNear(const HepLorentzVector & w,
00058 double epsilon) const {
00059 double limit = std::fabs(pp.dot(w.pp));
00060 limit += .25*((ee+w.ee)*(ee+w.ee));
00061 limit *= epsilon*epsilon;
00062 double delta = (pp - w.pp).mag2();
00063 delta += (ee-w.ee)*(ee-w.ee);
00064 return (delta <= limit );
00065 }
00066
00067 double HepLorentzVector::howNear(const HepLorentzVector & w) const {
00068 double wdw = std::fabs(pp.dot(w.pp)) + .25*((ee+w.ee)*(ee+w.ee));
00069 double delta = (pp - w.pp).mag2() + (ee-w.ee)*(ee-w.ee);
00070 if ( (wdw > 0) && (delta < wdw) ) {
00071 return std::sqrt (delta/wdw);
00072 } else if ( (wdw == 0) && (delta == 0) ) {
00073 return 0;
00074 } else {
00075 return 1;
00076 }
00077 }
00078
00079
00080
00081
00082
00083
00084 bool HepLorentzVector::isNearCM
00085 (const HepLorentzVector & w, double epsilon) const {
00086
00087 double tTotal = (ee + w.ee);
00088 Hep3Vector vTotal (pp + w.pp);
00089 double vTotal2 = vTotal.mag2();
00090
00091 if ( vTotal2 >= tTotal*tTotal ) {
00092
00093
00094
00095
00096 return (*this == w);
00097 }
00098
00099 if ( vTotal2 == 0 ) {
00100 return (isNear(w, epsilon));
00101 }
00102
00103
00104
00105 double tRecip = 1./tTotal;
00106 Hep3Vector bboost ( vTotal * (-tRecip) );
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116 double b2 = vTotal2*tRecip*tRecip;
00117
00118 register double ggamma = std::sqrt(1./(1.-b2));
00119 register double boostDotV1 = bboost.dot(pp);
00120 register double gm1_b2 = (ggamma-1)/b2;
00121
00122 HepLorentzVector w1 ( pp + ((gm1_b2)*boostDotV1+ggamma*ee) * bboost,
00123 ggamma * (ee + boostDotV1) );
00124
00125 register double boostDotV2 = bboost.dot(w.pp);
00126 HepLorentzVector w2 ( w.pp + ((gm1_b2)*boostDotV2+ggamma*w.ee) * bboost,
00127 ggamma * (w.ee + boostDotV2) );
00128
00129 return (w1.isNear(w2, epsilon));
00130
00131 }
00132
00133 double HepLorentzVector::howNearCM(const HepLorentzVector & w) const {
00134
00135 double tTotal = (ee + w.ee);
00136 Hep3Vector vTotal (pp + w.pp);
00137 double vTotal2 = vTotal.mag2();
00138
00139 if ( vTotal2 >= tTotal*tTotal ) {
00140
00141
00142
00143
00144 if (*this == w) {
00145 return 0;
00146 } else {
00147 return 1;
00148 }
00149 }
00150
00151 if ( vTotal2 == 0 ) {
00152 return (howNear(w));
00153 }
00154
00155
00156
00157 double tRecip = 1./tTotal;
00158 Hep3Vector bboost ( vTotal * (-tRecip) );
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168 double b2 = vTotal2*tRecip*tRecip;
00169
00170
00171
00172
00173 register double ggamma = std::sqrt(1./(1.-b2));
00174 register double boostDotV1 = bboost.dot(pp);
00175 register double gm1_b2 = (ggamma-1)/b2;
00176
00177 HepLorentzVector w1 ( pp + ((gm1_b2)*boostDotV1+ggamma*ee) * bboost,
00178 ggamma * (ee + boostDotV1) );
00179
00180 register double boostDotV2 = bboost.dot(w.pp);
00181 HepLorentzVector w2 ( w.pp + ((gm1_b2)*boostDotV2+ggamma*w.ee) * bboost,
00182 ggamma * (w.ee + boostDotV2) );
00183
00184 return (w1.howNear(w2));
00185
00186 }
00187
00188
00189
00190
00191
00192
00193
00194
00195 double HepLorentzVector::deltaR ( const HepLorentzVector & w ) const {
00196
00197 double a = eta() - w.eta();
00198 double b = pp.deltaPhi(w.getV());
00199
00200 return std::sqrt ( a*a + b*b );
00201
00202 }
00203
00204
00205
00206
00207
00208 bool HepLorentzVector::isParallel (const HepLorentzVector & w, double epsilon) const {
00209 double norm = euclideanNorm();
00210 double wnorm = w.euclideanNorm();
00211 if ( norm == 0 ) {
00212 if ( wnorm == 0 ) {
00213 return true;
00214 } else {
00215 return false;
00216 }
00217 }
00218 if ( wnorm == 0 ) {
00219 return false;
00220 }
00221 HepLorentzVector w1 = *this / norm;
00222 HepLorentzVector w2 = w / wnorm;
00223 return ( (w1-w2).euclideanNorm2() <= epsilon*epsilon );
00224 }
00225
00226
00227 double HepLorentzVector::howParallel (const HepLorentzVector & w) const {
00228
00229 double norm = euclideanNorm();
00230 double wnorm = w.euclideanNorm();
00231 if ( norm == 0 ) {
00232 if ( wnorm == 0 ) {
00233 return 0;
00234 } else {
00235 return 1;
00236 }
00237 }
00238 if ( wnorm == 0 ) {
00239 return 1;
00240 }
00241
00242 HepLorentzVector w1 = *this / norm;
00243 HepLorentzVector w2 = w / wnorm;
00244 double x1 = (w1-w2).euclideanNorm();
00245 return (x1 < 1) ? x1 : 1;
00246
00247 }
00248
00249 double HepLorentzVector::howLightlike() const {
00250 double m1 = std::fabs(restMass2());
00251 double twoT2 = 2*ee*ee;
00252 if (m1 < twoT2) {
00253 return m1/twoT2;
00254 } else {
00255 return 1;
00256 }
00257 }
00258
00259 }