00001 // @(#)root/physics:$Id: TLorentzVector.cxx 35182 2010-09-07 09:05:41Z brun $ 00002 // Author: Pasha Murat , Peter Malzacher 12/02/99 00003 // Oct 8 1999: changed Warning to Error and 00004 // return fX in Double_t & operator() 00005 // Oct 20 1999: dito in Double_t operator() 00006 // Jan 25 2000: implemented as (fP,fE) instead of (fX,fY,fZ,fE) 00007 00008 //______________________________________________________________________________ 00009 //*-*-*-*-*-*-*-*-*-*-*-*The Physics Vector package *-*-*-*-*-*-*-*-*-*-*-* 00010 //*-* ========================== * 00011 //*-* The Physics Vector package consists of five classes: * 00012 //*-* - TVector2 * 00013 //*-* - TVector3 * 00014 //*-* - TRotation * 00015 //*-* - TLorentzVector * 00016 //*-* - TLorentzRotation * 00017 //*-* It is a combination of CLHEPs Vector package written by * 00018 //*-* Leif Lonnblad, Andreas Nilsson and Evgueni Tcherniaev * 00019 //*-* and a ROOT package written by Pasha Murat. * 00020 //*-* for CLHEP see: http://wwwinfo.cern.ch/asd/lhc++/clhep/ * 00021 //*-* Adaption to ROOT by Peter Malzacher * 00022 //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* 00023 //BEGIN_HTML <!-- 00024 /* --> 00025 <H2>TLorentzVector</H2> 00026 <TT>TLorentzVector</TT> is a general four-vector class, which can be used 00027 either for the description of position and time (x,y,z,t) or momentum and 00028 energy (px,py,pz,E). 00029 <BR> 00030 <H3> 00031 Declaration</H3> 00032 <TT>TLorentzVector</TT> has been implemented as a set a <TT>TVector3</TT> and a <TT>Double_t</TT> variable. 00033 By default all components are initialized by zero. 00034 00035 <P><TT> TLorentzVector v1; // initialized 00036 by (0., 0., 0., 0.)</TT> 00037 <BR><TT> TLorentzVector v2(1., 1., 1., 1.);</TT> 00038 <BR><TT> TLorentzVector v3(v1);</TT> 00039 <BR><TT> TLorentzVector v4(TVector3(1., 2., 3.),4.);</TT> 00040 00041 <P>For backward compatibility there are two constructors from an <TT>Double_t</TT> 00042 and <TT>Float_t</TT> C array. 00043 <BR> 00044 00045 <H3> 00046 Access to the components</H3> 00047 There are two sets of access functions to the components of a<TT> LorentzVector</TT>: 00048 <TT>X(</TT>), <TT>Y()</TT>, <TT>Z()</TT>, <TT>T()</TT> and <TT>Px()</TT>,<TT> 00049 Py()</TT>, <TT>Pz()</TT> and <TT>E()</TT>. Both sets return the same values 00050 but the first set is more relevant for use where <TT>TLorentzVector</TT> 00051 describes a combination of position and time and the second set is more 00052 relevant where <TT>TLorentzVector</TT> describes momentum and energy: 00053 00054 <P><TT> Double_t xx =v.X();</TT> 00055 <BR><TT> ...</TT> 00056 <BR><TT> Double_t tt = v.T();</TT> 00057 00058 <P><TT> Double_t px = v.Px();</TT> 00059 <BR><TT> ...</TT> 00060 <BR><TT> Double_t ee = v.E();</TT> 00061 00062 <P>The components of TLorentzVector can also accessed by index: 00063 00064 <P><TT> xx = v(0); or 00065 xx = v[0];</TT> 00066 <BR><TT> yy = v(1); 00067 yy = v[1];</TT> 00068 <BR><TT> zz = v(2); 00069 zz = v[2];</TT> 00070 <BR><TT> tt = v(3); 00071 tt = v[3];</TT> 00072 00073 <P>You can use the <TT>Vect()</TT> member function to get the vector component 00074 of <TT>TLorentzVector</TT>: 00075 00076 <P> <TT>TVector3 p = v.Vect();</TT> 00077 00078 <P>For setting components also two sets of member functions can be used: 00079 SetX(),.., SetPx(),..: 00080 <BR> 00081 <BR><TT> v.SetX(1.); or 00082 v.SetPx(1.);</TT> 00083 <BR><TT> ... 00084 ...</TT> 00085 <BR><TT> v.SetT(1.); 00086 v.SetE(1.);</TT> 00087 00088 <P>To set more the one component by one call you can use the <TT>SetVect()</TT> 00089 function for the <TT>TVector3</TT> part or <TT>SetXYZT()</TT>, <TT>SetPxPyPzE()</TT>. For convenience there is 00090 00091 also a <TT>SetXYZM()</TT>: 00092 00093 <P><TT> v.SetVect(TVector3(1,2,3));</TT> 00094 <BR><TT> v.SetXYZT(x,y,z,t);</TT> 00095 <BR><TT> v.SetPxPyPzE(px,py,pz,e);</TT> 00096 <BR><TT> v.SetXYZM(x,y,z,m); // -> 00097 v=(x,y,z,e=Sqrt(x*x+y*y+z*z+m*m))</TT> 00098 <H3> 00099 Vector components in noncartesian coordinate systems</H3> 00100 There are a couple of memberfunctions to get and set the <TT>TVector3</TT> 00101 part of the parameters in 00102 <BR><B>sherical</B> coordinate systems: 00103 00104 <P><TT> Double_t m, theta, cost, phi, pp, pp2, ppv2, pp2v2;</TT> 00105 <BR><TT> m = v.Rho();</TT> 00106 <BR><TT> t = v.Theta();</TT> 00107 <BR><TT> cost = v.CosTheta();</TT> 00108 <BR><TT> phi = v.Phi();</TT> 00109 00110 <P><TT> v.SetRho(10.);</TT> 00111 <BR><TT> v.SetTheta(TMath::Pi()*.3);</TT> 00112 <BR><TT> v.SetPhi(TMath::Pi());</TT> 00113 00114 <P>or get infoormation about the r-coordinate in <B>cylindrical</B> systems: 00115 00116 <P><TT> Double_t pp, pp2, ppv2, pp2v2;</TT> 00117 <BR><TT> pp = v.Perp(); // get transvers component</TT> 00118 <BR><TT> pp2 = v.Perp2(); // get transverse component squared</TT> 00119 <BR><TT> ppv2 = v.Perp(v1); // get 00120 transvers component with</TT> 00121 <BR><TT> // respect to another vector</TT> 00122 <BR><TT> pp2v2 = v.Perp(v1);</TT> 00123 00124 <P>for convenience there are two more set functions <TT>SetPtEtaPhiE(pt,eta,phi,e); 00125 and SetPtEtaPhiM(pt,eta,phi,m);</TT> 00126 <H3> 00127 Arithmetic and comparison operators</H3> 00128 The <TT>TLorentzVecto</TT>r class provides operators to add, subtract or 00129 compare four-vectors: 00130 00131 <P><TT> v3 = -v1;</TT> 00132 <BR><TT> v1 = v2+v3;</TT> 00133 <BR><TT> v1+= v3;</TT> 00134 <BR><TT> v1 = v2 + v3;</TT> 00135 <BR><TT> v1-= v3;</TT> 00136 00137 <P><TT> if (v1 == v2) {...}</TT> 00138 <BR><TT> if(v1 != v3) {...}</TT> 00139 <H3> 00140 Magnitude/Invariant mass, beta, gamma, scalar product</H3> 00141 The scalar product of two four-vectors is calculated with the (-,-,-,+) 00142 metric, 00143 <BR> i.e. <TT><B>s = v1*v2 = </B>t1*t2-x1*x2-y1*y2-z1*z2</TT> 00144 <BR>The magnitude squared <B><TT>mag2</TT></B> of a four-vector is therfore: 00145 <BR> <TT>mag2<B> 00146 = v*v = </B>t*t-x*x-y*y-z*z</TT> 00147 <BR>It <TT>mag2</TT> is negative <TT>mag = -Sqrt(-mag*mag)</TT>. The member 00148 functions are: 00149 00150 <P><TT> Double_t s, s2;</TT> 00151 <BR><TT> s = v1.Dot(v2); // scalar 00152 product</TT> 00153 <BR><TT> s = v1*v2; // scalar product</TT> 00154 <BR><TT> s2 = v.Mag2(); or s2 = v.M2();</TT> 00155 <BR><TT> s = v.Mag(); 00156 s = v.M();</TT> 00157 00158 <P>Since in case of momentum and energy the magnitude has the meaning of 00159 invariant mass <TT>TLorentzVector</TT> provides the more meaningful aliases 00160 <TT>M2()</TT> and <TT>M()</TT>; 00161 <P>The member functions <TT>Beta()</TT> and <TT>Gamma()</TT> returns 00162 <TT>beta</TT> and <tt>gamma = 1/Sqrt(1-beta*beta)</tt>. 00163 <H3> 00164 Lorentz boost</H3> 00165 A boost in a general direction can be parameterized with three parameters 00166 which can be taken as the components of a three vector <TT><B>b </B>= (bx,by,bz)</TT>. 00167 With 00168 <BR><TT><B> x</B> = (x,y,z) and gamma = 1/Sqrt(1-beta*beta)</TT> (beta being the module of vector b), 00169 an arbitary active Lorentz boost transformation (from the rod frame 00170 to the original frame) can be written as: 00171 <BR> <TT><B>x</B> 00172 = <B>x'</B> + (gamma-1)/(beta*beta) * (<B>b</B>*<B>x'</B>) * <B>b</B> + 00173 gamma * t' *<B> b</B></TT> 00174 <BR><B> </B><TT>t 00175 = gamma (t'+ <B>b</B>*<B>x'</B>).</TT> 00176 00177 <P>The member function <TT>Boost()</TT> performs a boost transformation 00178 from the rod frame to the original frame. <TT>BoostVector()</TT> returns 00179 a <TT>TVector3</TT> of the spatial components divided by the time component: 00180 00181 <P><TT> TVector3 b;</TT> 00182 <BR><TT> v.Boost(bx,by,bz);</TT> 00183 <BR><TT> v.Boost(b);</TT> 00184 <BR><TT> b = v.BoostVector(); // b=(x/t,y/t,z/t)</TT> 00185 <H3> 00186 Rotations</H3> 00187 There are four sets of functions to rotate the <TT>TVector3</TT> component 00188 of a <TT>TLorentzVector</TT>: 00189 <H5> 00190 rotation around axes</H5> 00191 <TT> v.RotateX(TMath::Pi()/2.);</TT> 00192 <BR><TT> v.RotateY(.5);</TT> 00193 <BR><TT> v.RotateZ(.99);</TT> 00194 <H5> 00195 rotation around an arbitary axis</H5> 00196 <TT> v.Rotate(TMath::Pi()/4., v1); // rotation around v1</TT> 00197 <H5> 00198 transformation from rotated frame</H5> 00199 <TT> v.RotateUz(direction); // direction must be a unit TVector3</TT> 00200 <H5> 00201 by TRotation (see TRotation)</H5> 00202 <TT> TRotation r;</TT> 00203 <BR><TT> v.Transform(r); or 00204 v *= r; // Attention v=M*v</TT> 00205 <H3> 00206 Misc</H3> 00207 00208 <H5> 00209 Angle between two vectors</H5> 00210 <TT> Double_t a = v1.Angle(v2.Vect()); // get angle between v1 and 00211 v2</TT> 00212 <H5> 00213 Light-cone components</H5> 00214 Member functions <TT>Plus()</TT> and <TT>Minus(</TT>) return the positive 00215 and negative light-cone components:<TT></TT> 00216 00217 <P><TT> Double_t pcone = v.Plus();</TT> 00218 <BR><TT> Double_t mcone = v.Minus();</TT> 00219 <P>CAVEAT: The values returned are T{+,-}Z. It is known that some authors 00220 find it easier to define these components as (T{+,-}Z)/sqrt(2). Thus 00221 check what definition is used in the physics you're working in and adapt 00222 your code accordingly. 00223 00224 <H5> 00225 Transformation by TLorentzRotation</H5> 00226 A general Lorentz transformation see class <TT>TLorentzRotation</TT> can 00227 be used by the <TT>Transform()</TT> member function, the <TT>*=</TT> or 00228 <TT>*</TT> operator of the <TT>TLorentzRotation</TT> class:<TT></TT> 00229 00230 <P><TT> TLorentzRotation l;</TT> 00231 <BR><TT> v.Transform(l);</TT> 00232 <BR><TT> v = l*v; or 00233 v *= l; // Attention v = l*v</TT> 00234 00235 <P> 00236 <!--*/ 00237 // -->END_HTML 00238 00239 #include "TClass.h" 00240 #include "TError.h" 00241 #include "TLorentzVector.h" 00242 #include "TLorentzRotation.h" 00243 00244 ClassImp(TLorentzVector) 00245 00246 TLorentzVector::TLorentzVector(Double_t x, Double_t y, Double_t z, Double_t t) 00247 : fP(x,y,z), fE(t) {} 00248 00249 TLorentzVector::TLorentzVector(const Double_t * x0) 00250 : fP(x0), fE(x0[3]) {} 00251 00252 TLorentzVector::TLorentzVector(const Float_t * x0) 00253 : fP(x0), fE(x0[3]) {} 00254 00255 TLorentzVector::TLorentzVector(const TVector3 & p, Double_t e) 00256 : fP(p), fE(e) {} 00257 00258 TLorentzVector::TLorentzVector(const TLorentzVector & p) : TObject(p) 00259 , fP(p.Vect()), fE(p.T()) {} 00260 00261 TLorentzVector::~TLorentzVector() {} 00262 00263 Double_t TLorentzVector::operator () (int i) const 00264 { 00265 //dereferencing operator const 00266 switch(i) { 00267 case kX: 00268 case kY: 00269 case kZ: 00270 return fP(i); 00271 case kT: 00272 return fE; 00273 default: 00274 Error("operator()()", "bad index (%d) returning 0",i); 00275 } 00276 return 0.; 00277 } 00278 00279 Double_t & TLorentzVector::operator () (int i) 00280 { 00281 //dereferencing operator 00282 switch(i) { 00283 case kX: 00284 case kY: 00285 case kZ: 00286 return fP(i); 00287 case kT: 00288 return fE; 00289 default: 00290 Error("operator()()", "bad index (%d) returning &fE",i); 00291 } 00292 return fE; 00293 } 00294 00295 void TLorentzVector::Boost(Double_t bx, Double_t by, Double_t bz) 00296 { 00297 //Boost this Lorentz vector 00298 Double_t b2 = bx*bx + by*by + bz*bz; 00299 register Double_t gamma = 1.0 / TMath::Sqrt(1.0 - b2); 00300 register Double_t bp = bx*X() + by*Y() + bz*Z(); 00301 register Double_t gamma2 = b2 > 0 ? (gamma - 1.0)/b2 : 0.0; 00302 00303 SetX(X() + gamma2*bp*bx + gamma*bx*T()); 00304 SetY(Y() + gamma2*bp*by + gamma*by*T()); 00305 SetZ(Z() + gamma2*bp*bz + gamma*bz*T()); 00306 SetT(gamma*(T() + bp)); 00307 } 00308 00309 Double_t TLorentzVector::Rapidity() const 00310 { 00311 //return rapidity 00312 return 0.5*log( (E()+Pz()) / (E()-Pz()) ); 00313 } 00314 00315 TLorentzVector &TLorentzVector::operator *= (const TLorentzRotation & m) 00316 { 00317 //multiply this Lorentzvector by m 00318 return *this = m.VectorMultiplication(*this); 00319 } 00320 00321 TLorentzVector &TLorentzVector::Transform(const TLorentzRotation & m) 00322 { 00323 //Transform this Lorentzvector 00324 return *this = m.VectorMultiplication(*this); 00325 } 00326 00327 void TLorentzVector::Streamer(TBuffer &R__b) 00328 { 00329 // Stream an object of class TLorentzVector. 00330 Double_t x, y, z; 00331 UInt_t R__s, R__c; 00332 if (R__b.IsReading()) { 00333 Version_t R__v = R__b.ReadVersion(&R__s, &R__c); 00334 if (R__v > 3) { 00335 R__b.ReadClassBuffer(TLorentzVector::Class(), this, R__v, R__s, R__c); 00336 return; 00337 } 00338 //====process old versions before automatic schema evolution 00339 if (R__v != 2) TObject::Streamer(R__b); 00340 R__b >> x; 00341 R__b >> y; 00342 R__b >> z; 00343 fP.SetXYZ(x,y,z); 00344 R__b >> fE; 00345 R__b.CheckByteCount(R__s, R__c, TLorentzVector::IsA()); 00346 } else { 00347 R__b.WriteClassBuffer(TLorentzVector::Class(),this); 00348 } 00349 } 00350 00351 00352 //______________________________________________________________________________ 00353 void TLorentzVector::Print(Option_t *) const 00354 { 00355 // Print the TLorentz vector components as (x,y,z,t) and (P,eta,phi,E) representations 00356 Printf("(x,y,z,t)=(%f,%f,%f,%f) (P,eta,phi,E)=(%f,%f,%f,%f)", 00357 fP.x(),fP.y(),fP.z(),fE, 00358 P(),Eta(),Phi(),fE); 00359 }