#ifndef INCLUDE_TMBVECTOR3 #define INCLUDE_TMBVECTOR3 //______________________________________________________________________________ //*-*-*-*-*-*-*-*-*-*-*-*The Physics Vector package *-*-*-*-*-*-*-*-*-*-*-* //*-* ========================== * //*-* The Physics Vector package consists of five classes: * //*-* - TVector2 * //*-* - TVector3 * //*-* - TRotation * //*-* - TLorentzVector * //*-* - TLorentzRotation * //*-* It is a combination of CLHEPs Vector package written by * //*-* Leif Lonnblad, Andreas Nilsson and Evgueni Tcherniaev * //*-* and a ROOT package written by Pasha Murat. * //*-* for CLHEP see: http://wwwinfo.cern.ch/asd/lhc++/clhep/ * //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* //BEGIN_HTML
TVector3 v1; //
v1 = (0,0,0)
TVector3 v2(1); // v2 = (1,0,0)
TVector3 v3(1,2,3); // v3 = (1,2,3)
TVector3 v4(v2); // v4 = v2
It is also possible (but not recommended) to initialize a TVector3 with a Double_t or Float_t C array.
You can get the basic components either by name or by index using operator():
xx = v1.X(); or xx =
v1(0);
yy = v1.Y();
yy = v1(1);
zz = v1.Z();
zz = v1(2);
The memberfunctions SetX(), SetY(), SetZ() and SetXYZ() allow to set the components:
v1.SetX(1.); v1.SetY(2.); v1.SetZ(3.);
v1.SetXYZ(1.,2.,3.);
Double_t m = v.Mag3(); // get magnitude
(=rho=Sqrt(x*x+y*y+z*z)))
Double_t m2 = v.Mag32(); // get magnitude squared
Double_t t = v.Theta(); // get polar angle
Double_t ct = v.CosTheta();// get cos of theta
Double_t p = v.Phi(); // get azimuth angle
Double_t pp = v.Perp(); // get transverse component
Double_t pp2= v.Perp2(); // get transvers component
squared
It is also possible to get the transverse component with respect to another vector:
Double_t ppv1 = v.Perp(v1);
Double_t pp2v1 = v.Perp2(v1);
The pseudorapiditiy ( eta=-ln (tan (phi/2)) ) can be get by Eta()
or PseudoRapidity():
Double_t eta = v.PseudoRapidity();
There are set functions to change one of the noncartesian coordinates:
v.SetTheta(.5); // keeping rho and phi
v.SetPhi(.8); // keeping rho and theta
v.SetMag3(10.); // keeping theta and phi
v.SetPerp(3.); // keeping z and phi
v3 = -v1;
v1 = v2+v3;
v1 += v3;
v1 = v1 - v3
v1 -= v3;
v1 *= 10;
v1 = 5*v2;
if(v1==v2) {...}
if(v1!=v2) {...}
TRotation m;
...
v1.transform(m);
v1 = m*v1;
v1 *= m; // Attention v1 = m*v1
transforms v1 from the rotated frame (z' parallel to direction, x' in
the theta plane and y' in the xy plane as well as perpendicular to the
theta plane) to the (x,y,z) frame.'
END_HTML
//
// See copyright statements at end of file.
#include "TMath.h"
#include "TError.h"
#include "TVector2.h"
#include "TMatrix.h"
#include "TRotation.h"
class TMBVector3 : public TObject {
public:
TMBVector3(Double_t x = 0.0, Double_t y = 0.0, Double_t z = 0.0):
fX(x), fY(y), fZ(z) { /// The constructor.
}
TMBVector3(const Double_t *a): fX(a[0]), fY(a[1]), fZ(a[2])
{ /// Constructor from an array
}
TMBVector3(const Float_t *a): fX(a[0]), fY(a[1]), fZ(a[2])
{ /// Constructor from an array
}
TMBVector3(const TMBVector3 &p): TObject(p),
fX(p.fX), fY(p.fY), fZ(p.fZ)
{ /// Copy constructor.
}
TMBVector3(const TVector3& v): fX(v.X()), fY(v.Y()), fZ(v.Z())
{
/// conversion constructor from TVector3
}
virtual ~TMBVector3() {}
// Destructor
Double_t operator () (int i) const {
/// Get components by index.
switch(i) {
case 0: return fX;
case 1: return fY;
case 2: return fZ;
default: Error("operator()(i)", "bad index (%d) returning &fX",i);
}
return 0.;
}
inline Double_t operator [] (int i) const { /// Get components by index.
return operator()(i); }
Double_t & operator () (int i) {
/// Get components by index.
switch(i) {
case 0: return fX;
case 1: return fY;
case 2: return fZ;
default: Error("operator()(i)", "bad index (%d) returning &fX",i);
}
return fX;
}
inline Double_t & operator [] (int i) { /// Set components by index.
return operator()(i); }
Double_t x() const { /// get X component
return fX; }
Double_t X() const { /// get X component
return x(); }
Double_t Px() const { /// get X component
return x(); }
Double_t y() const { /// get Y component
return fY; }
Double_t Y() const { /// get Y component
return y(); }
Double_t Py() const { /// get X component
return y(); }
Double_t z() const { /// get Y component
return fZ; }
Double_t Z() const { /// get Z component
return z(); }
Double_t Pz() const { /// get X component
return z(); }
void SetX(Double_t a) { /// set X component, see TVector3::SetX()
fX=a; }
void SetPx(Double_t a) { /// set X component, see TVector3::SetX()
SetX(a); }
void SetY(Double_t a) { /// set X component, see TVector3::SetY()
fY=a; }
void SetPy(Double_t a) { /// set X component, see TVector3::SetX()
SetY(a); }
void SetZ(Double_t a) { /// set X component, see TVector3::SetZ()
fZ=a; }
void SetPz(Double_t a) { /// set X component, see TVector3::SetX()
SetZ(a); }
inline void SetXYZ(Double_t x, Double_t y, Double_t z) {
/// set all members using Cartesian coordinates
fX=x; fY=y; fZ=z; }
void SetPtEtaPhi(Double_t pt, Double_t eta, Double_t phi) {
/// set all members using |momentum|, eta, and phi
Double_t apt = TMath::Abs(pt);
Double_t d=TMath::Tan(2.0*TMath::ATan(TMath::Exp(-eta)));
if (d==0.)
Warning("SetPtEtaPhi", "tan(2*atan(exp(-eta)))==0.! Setting z=0.");
SetXYZ(apt*TMath::Cos(phi), apt*TMath::Sin(phi), d==0.?0.:apt/d );
}
void SetPtThetaPhi(Double_t pt, Double_t theta, Double_t phi) {
/// set all members using |momentum|, theta, and phi
Double_t apt = TMath::Abs(pt);
fX = apt * TMath::Cos(phi);
fY = apt * TMath::Sin(phi);
Double_t tanTheta = TMath::Tan(theta);
if (tanTheta==0.)
Warning("SetPtThetaPhi", "tan(theta)==0.! Setting z=0.");
fZ = tanTheta ? apt / tanTheta : 0;
}
inline void GetXYZ(Double_t *carray) const {
/// Getters into an arry (no checking of array bounds!)
carray[0]=fX; carray[1]=fY; carray[2]=fZ; }
inline void GetXYZ(Float_t *carray) const {
/// Getters into an arry (no checking of array bounds!)
carray[0]=fX; carray[1]=fY; carray[2]=fZ; }
Double_t Eta() const { /// get eta (aka pseudo-rapidity)
return PseudoRapidity(); }
Double_t Theta() const { /// get theta
return fX==.0 && fY==.0 && fZ==.0 ? .0 : TMath::ATan2(Perp(),fZ); }
Double_t CosTheta() const { /// get cos(Theta)
Double_t ptot = Mag3();
return ptot == 0.0 ? 1.0 : fZ/ptot;
}
Double_t Phi() const { /// get phi from 0..2Pi
return (fX==.0 && fY==.0 ? .0 : TVector2::Phi_0_2pi(TMath::ATan2(fY,fX))); }
Double_t Rho() const { /// get magnitude of momentum ("radius" in spherical coords)
return Mag3(); }
virtual Double_t Mag32() const {
/// The magnitude squared (rho^2 in spherical coordinate system).
return fX*fX + fY*fY + fZ*fZ; }
virtual Double_t Mag3() const {
/// The magnitude (rho in spherical coordinate system).
return TMath::Sqrt(Mag32()); }
Double_t P() const { /// get |P|
return Mag3(); }
Double_t DeltaPhi(const TMBVector3 &v) const {
/// Get difference in phi, between -pi and pi
return TVector2::Phi_mpi_pi(Phi()-v.Phi()); }
Double_t DeltaR(const TMBVector3 &v) const {
/// Get "distance" dR of v to *this, dR=sqrt(dPhi*dPhi+dEta*dEta)
Double_t deta = Eta()-v.Eta();
Double_t dphi = TVector2::Phi_mpi_pi(Phi()-v.Phi());
return TMath::Sqrt( deta*deta+dphi*dphi );
}
Double_t DrEtaPhi(const TMBVector3 &lv) const {
/// Get "distance" dR of lv to *this, dR=sqrt(dPhi*dPhi+dEta*dEta)
return DeltaR(lv); }
// TVector2 EtaPhiVector() { /// return TVector2(eta,phi)
// return TVector2 (Eta(),Phi()); }
Double_t Angle(const TVector3 & q) const { /// Angle wrt. another vector.
Double_t ptot2 = Mag32()*q.Mag2();
if(ptot2 <= 0.) return .0;
Double_t arg = Dot(q)/TMath::Sqrt(ptot2);
if(arg > 1.0) arg = 1.0;
if(arg < -1.0) arg = -1.0;
return TMath::ACos(arg);
}
void SetPhi(Double_t ph) { /// set phi keeping mag and theta constant
Double_t xy = Perp();
SetX(xy*TMath::Cos(ph));
SetY(xy*TMath::Sin(ph));
}
inline void SetTheta(Double_t th) { /// Set theta keeping mag and phi constant
Double_t ma = Mag3();
Double_t ph = Phi();
SetX(ma*TMath::Sin(th)*TMath::Cos(ph));
SetY(ma*TMath::Sin(th)*TMath::Sin(ph));
SetZ(ma*TMath::Cos(th));
}
inline void SetMag3(Double_t ma) { /// Set magnitude keeping theta and phi constant
Double_t factor = Mag3();
if (factor == 0) {
Warning("SetMag3","zero vector can't be stretched");
}else{
factor = ma/factor;
SetX(fX*factor);
SetY(fY*factor);
SetZ(fZ*factor);
}
}
inline Double_t Perp2() const {
/// The transverse component squared (R^2 in cylindrical coordinate system)
return fX*fX + fY*fY; }
inline Double_t Pt() const {
/// transverse component; projection of 3 vector onto XY plane (R in cylindrical coords)
return Perp(); }
inline Double_t Perp() const {
/// The transverse component (R in cylindrical coordinate system)
return TMath::Sqrt(Perp2()); }
inline void SetPerp(Double_t r) {
/// Set the transverse component keeping phi and z constant.
Double_t p = Perp();
if (p != 0.0) {
fX *= r/p;
fY *= r/p;
}
}
inline Double_t Perp2(const TMBVector3 &p) const {
/// The transverse component w.r.t. given axis squared.
Double_t tot = p.Mag32();
Double_t ss = Dot(p);
return tot > 0.0 ? Mag32()-ss*ss/tot : Mag32();
}
inline Double_t Pt(const TMBVector3 &p) const {
/// The transverse component w.r.t. given axis.
return Perp(p); }
inline Double_t Perp(const TMBVector3 &p) const {
/// The transverse component w.r.t. given axis.
return TMath::Sqrt(Perp2(p)); }
inline void SetMag3ThetaPhi(Double_t mag, Double_t theta, Double_t phi) {
/// Set all components as magnitude, theta, and phi
Double_t amag = TMath::Abs(mag);
fX = amag * TMath::Sin(theta) * TMath::Cos(phi);
fY = amag * TMath::Sin(theta) * TMath::Sin(phi);
fZ = amag * TMath::Cos(theta);
}
inline TMBVector3 & operator = (const TMBVector3 &p) {
/// Assignment.
fX = p.fX; fY = p.fY; fZ = p.fZ;
return *this;
}
inline Bool_t operator == (const TMBVector3 &v) const {
/// Equality operator. Equal if all components X,Y,Z are equal.
return (v.fX==fX && v.fY==fY && v.fZ==fZ); }
inline Bool_t operator != (const TMBVector3 &v) const {
/// Inequality operator. Not equal if any of X,Y,Z are not equal.
return (v.fX!=fX || v.fY!=fY || v.fZ!=fZ); }
/// Equivalence operator. True if all components are
/// equivalent within machine precision.
Bool_t is_equal (const TMBVector3 &v) const;
inline TMBVector3 & operator += (const TMBVector3 &p) { /// Addition.
fX += p.fX; fY += p.fY; fZ += p.fZ;
return *this;
}
inline TMBVector3 & operator -= (const TMBVector3 &p) { /// Subtraction.
fX -= p.fX; fY -= p.fY; fZ -= p.fZ;
return *this;
}
inline TMBVector3 operator - () const { /// Unary minus.
return TMBVector3(-fX, -fY, -fZ); }
inline TMBVector3 & operator *= (Double_t a) {
/// Scaling with real numbers.
fX *= a; fY *= a; fZ *= a;
return *this;
}
inline TMBVector3 Unit() const { /// Unit vector parallel to this.
Double_t tot = Mag32();
TMBVector3 p(fX,fY,fZ);
return tot>.0 ? p *= (1.0/TMath::Sqrt(tot)) : p;
}
inline TMBVector3 Orthogonal() const { /// Vector orthogonal to this.
Double_t x = fX < 0.0 ? -fX : fX;
Double_t y = fY < 0.0 ? -fY : fY;
Double_t z = fZ < 0.0 ? -fZ : fZ;
if (x < y)
return x