MultivariateAnalysis/interface/TMBVector3.h

#ifndef INCLUDE_TMBVECTOR3
#define INCLUDE_TMBVECTOR3

#include "TMath.h"
#include "TError.h"
#include "TVector2.h"
#include "TMatrix.h"
#include "TRotation.h"

/*! \brief 
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/           

\ingroup reco

<H2>
TVector3</H2>
<TT>TVector3</TT> is a general three vector class, which can be used for
the description of different vectors in 3D.
<H3>
Declaration / Access to the components</H3>
<TT>TVector3</TT> has been implemented as a vector of three <TT>Double_t</TT>
variables, representing the cartesian coordinates. By default all components
are initialized to zero:

<P><TT>&nbsp; TVector3 v1;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; //
v1 = (0,0,0)</TT>
<BR><TT>&nbsp; TVector3 v2(1);&nbsp;&nbsp;&nbsp;&nbsp; // v2 = (1,0,0)</TT>
<BR><TT>&nbsp; TVector3 v3(1,2,3); // v3 = (1,2,3)</TT>
<BR><TT>&nbsp; TVector3 v4(v2);&nbsp;&nbsp;&nbsp; // v4 = v2</TT>

<P>It is also possible (but not recommended) to initialize a <TT>TVector3</TT>
with a <TT>Double_t</TT> or <TT>Float_t</TT> C array.

<P>You can get the basic components either by name or by index using <TT>operator()</TT>:

<P><TT>&nbsp; xx = v1.X();&nbsp;&nbsp;&nbsp; or&nbsp;&nbsp;&nbsp; xx =
v1(0);</TT>
<BR><TT>&nbsp; yy = v1.Y();&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
yy = v1(1);</TT>
<BR><TT>&nbsp; zz = v1.Z();&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
zz = v1(2);</TT>

<P>The memberfunctions <TT>SetX()</TT>, <TT>SetY()</TT>, <TT>SetZ()</TT>
and<TT> SetXYZ()</TT> allow to set the components:

<P><TT>&nbsp; v1.SetX(1.); v1.SetY(2.); v1.SetZ(3.);</TT>
<BR><TT>&nbsp; v1.SetXYZ(1.,2.,3.);</TT>
<BR>&nbsp;
<H3>
Noncartesian coordinates</H3>
To get information on the <TT>TVector3</TT> in spherical (rho,phi,theta)
or cylindrical (z,r,theta) coordinates, the
<BR>the member functions <TT>Mag3()</TT> (=magnitude=rho in spherical coordinates),
<TT>Mag32()</TT>, <TT>Theta()</TT>, <TT>CosTheta()</TT>, <TT>Phi()</TT>,
<TT>Perp()</TT> (the transverse component=r in cylindrical coordinates),
<TT>Perp2()</TT> can be used:

<P><TT>&nbsp; Double_t m&nbsp; = v.Mag3();&nbsp;&nbsp;&nbsp; // get magnitude
(=rho=Sqrt(x*x+y*y+z*z)))</TT>
<BR><TT>&nbsp; Double_t m2 = v.Mag32();&nbsp;&nbsp; // get magnitude squared</TT>
<BR><TT>&nbsp; Double_t t&nbsp; = v.Theta();&nbsp; // get polar angle</TT>
<BR><TT>&nbsp; Double_t ct = v.CosTheta();// get cos of theta</TT>
<BR><TT>&nbsp; Double_t p&nbsp; = v.Phi();&nbsp;&nbsp;&nbsp; // get azimuth angle</TT>
<BR><TT>&nbsp; Double_t pp = v.Perp();&nbsp;&nbsp; // get transverse component</TT>
<BR><TT>&nbsp; Double_t pp2= v.Perp2();&nbsp; // get transvers component
squared</TT>

<P>It is also possible to get the transverse component with respect to
another vector:

<P><TT>&nbsp; Double_t ppv1 = v.Perp(v1);</TT>
<BR><TT>&nbsp; Double_t pp2v1 = v.Perp2(v1);</TT>

<P>The pseudorapiditiy ( eta=-ln (tan (phi/2)) ) can be get by <TT>Eta()</TT>
or <TT>PseudoRapidity()</TT>:
<BR>&nbsp;
<BR><TT>&nbsp; Double_t eta = v.PseudoRapidity();</TT>

<P>There are set functions to change one of the noncartesian coordinates:

<P><TT>&nbsp; v.SetTheta(.5); // keeping rho and phi</TT>
<BR><TT>&nbsp; v.SetPhi(.8);&nbsp;&nbsp; // keeping rho and theta</TT>
<BR><TT>&nbsp; v.SetMag3(10.);&nbsp; // keeping theta and phi</TT>
<BR><TT>&nbsp; v.SetPerp(3.);&nbsp; // keeping z and phi</TT>
<BR>&nbsp;
<H3>
Arithmetic / Comparison</H3>
The <TT>TVector3</TT> class provides the operators to add, subtract, scale and compare
vectors:

<P><TT>&nbsp; v3&nbsp; = -v1;</TT>
<BR><TT>&nbsp; v1&nbsp; = v2+v3;</TT>
<BR><TT>&nbsp; v1 += v3;</TT>
<BR><TT>&nbsp; v1&nbsp; = v1 - v3</TT>
<BR><TT>&nbsp; v1 -= v3;</TT>
<BR><TT>&nbsp; v1 *= 10;</TT>
<BR><TT>&nbsp; v1&nbsp; = 5*v2;</TT>

<P><TT>&nbsp; if(v1==v2) {...}</TT>
<BR><TT>&nbsp; if(v1!=v2) {...}</TT>
<BR>&nbsp;
<H3>
Related Vectors</H3>
<TT>&nbsp; v2 = v1.Unit();&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; // get unit
vector parallel to v1</TT>
<BR><TT>&nbsp; v2 = v1.Orthogonal(); // get vector orthogonal to v1</TT>
<H3>
Scalar and vector products</H3>
<TT>&nbsp; s = v1.Dot(v2);&nbsp;&nbsp; // scalar product</TT>
<BR><TT>&nbsp; s = v1 * v2;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; // scalar product</TT>
<BR><TT>&nbsp; v = v1.Cross(v2); // vector product</TT>
<H3>
&nbsp;Angle between two vectors</H3>
<TT>&nbsp; Double_t a = v1.Angle(v2);</TT>
<H3>
Rotations</H3>

<H5>
Rotation around axes</H5>
<TT>&nbsp; v.RotateX(.5);</TT>
<BR><TT>&nbsp; v.RotateY(TMath::Pi());</TT>
<BR><TT>&nbsp; v.RotateZ(angle);</TT>
<H5>
Rotation around a vector</H5>
<TT>&nbsp; v1.Rotate(TMath::Pi()/4, v2); // rotation around v2</TT>
<H5>
Rotation by TRotation</H5>
<TT>TVector3</TT> objects can be rotated by objects of the <TT>TRotation</TT>
class using the <TT>Transform()</TT> member functions,
<BR>the <TT>operator *=</TT> or the <TT>operator *</TT> of the TRotation
class:

<P><TT>&nbsp; TRotation m;</TT>
<BR><TT>&nbsp; ...</TT>
<BR><TT>&nbsp; v1.transform(m);</TT>
<BR><TT>&nbsp; v1 = m*v1;</TT>
<BR><TT>&nbsp; v1 *= m; // Attention v1 = m*v1</TT>
<H5>
Transformation from rotated frame</H5>
<TT>&nbsp; TVector3 direction = v.Unit()</TT>
<BR><TT>&nbsp; v1.RotateUz(direction); // direction must be TVector3 of
unit length</TT>

<P>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.'

 See copyright statements at end of file.

*/

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<z ? TMBVector3(0,fZ,-fY) : TMBVector3(fY,-fX,0);
     return y<z ? TMBVector3(-fZ,0,fX) : TMBVector3(fY,-fX,0);
  }

  Double_t Dot(const TMBVector3 &p) const { /// Scalar product.
     return fX*p.fX + fY*p.fY + fZ*p.fZ; }

  inline TMBVector3 Cross(const TMBVector3 &p) const { /// Cross product
     return TMBVector3(fY*p.fZ-p.fY*fZ, fZ*p.fX-p.fZ*fX, fX*p.fY-p.fX*fY);
  }

  Double_t PseudoRapidity() const {
     /// Returns the pseudo-rapidity, i.e. -ln(tan(theta/2))
     double cosTheta = CosTheta();
     if (cosTheta*cosTheta < 1) 
        return (-0.5* TMath::Log( (1.0-cosTheta)/(1.0+cosTheta) ));
     Warning("PseudoRapidity","transvers momentum = 0! return +/- 10e10");
     if (fZ > 0) return (10e10);
     else        return (-10e10);
  }
  void RotateX(Double_t angle) {
     /// Rotates around the x-axis.
     Double_t s = TMath::Sin(angle);
     Double_t c = TMath::Cos(angle);
     Double_t yy = fY;
     fY = c*yy - s*fZ;
     fZ = s*yy + c*fZ;
  }

  void RotateY(Double_t angle) {
     /// Rotates around the y-axis.
     Double_t s = TMath::Sin(angle);
     Double_t c = TMath::Cos(angle);
     Double_t zz = fZ;
     fZ = c*zz - s*fX;
     fX = s*zz + c*fX;
  }

  void RotateZ(Double_t angle) {
     /// Rotates around the z-axis.
     Double_t s = TMath::Sin(angle);
     Double_t c = TMath::Cos(angle);
     Double_t xx = fX;
     fX = c*xx - s*fY;
     fY = s*xx + c*fY;
  }

  void RotateUz(const TMBVector3& NewUzVector) {
     /// Rotates reference frame from Uz to newUz (must be unit vector!)
     Double_t u1 = NewUzVector.fX;
     Double_t u2 = NewUzVector.fY;
     Double_t u3 = NewUzVector.fZ;
     Double_t up = u1*u1 + u2*u2;

     if (up>0.) {
        up = TMath::Sqrt(up);
        Double_t px = fX,  py = fY,  pz = fZ;
        fX = (u1*u3*px - u2*py + u1*up*pz)/up;
        fY = (u2*u3*px + u1*py + u2*up*pz)/up;
        fZ = (u3*u3*px -    px + u3*up*pz)/up;
     }
     else if (u3 < 0.) { fX = -fX; fZ = -fZ; }      // phi=0  teta=pi
     else {};
  }

  void Rotate(Double_t angle, const TMBVector3 &axis) {
     /// Rotates around the axis specified by another vector
     TRotation trans;
     trans.Rotate(angle, (TVector3)axis);
     operator*=(trans);
  }

  TMBVector3 & operator *= (const TRotation &m) {
     /// Transform with a rotation matrix
     return *this = m * (*this);
  }
  TMBVector3 & Transform(const TRotation &m) {
     /// Transform with a rotation matrix
     return *this = m * (*this);
  }

   TVector2 XYvector() const {
     /// TVector2 containing x and y components
     return TVector2(fX,fY); }

   operator TVector3() const { 
      /// cast to a TVector3
      return TVector3(X(),Y(),Z()); }

   void Print(Option_t* option="") const {
      /// print components to stdout
      Printf("%s %s (x,y,z)=(%f,%f,%f) (rho,theta,phi)=(%f,%f,%f)",
             GetName(),GetTitle(),X(),Y(),Z(),
             Mag3(),Theta()*TMath::RadToDeg(),Phi()*TMath::RadToDeg());
   }

   // Helper to test if two FP numbers are equivalent within machine precision.
   static bool is_equal (double x1, double x2);

private:
   Double32_t fX; ///< X component (left-right)
   Double32_t fY; ///< Y component (up-down)
   Double32_t fZ; ///< Z component (along the beam)

     //ClassDef(TMBVector3,1) // A 3D physics vector

};

inline
TMBVector3 operator + (const TMBVector3 &a, const TMBVector3 &b) {
   /// Addition of 3-vectors.
  return TVector3(a.X() + b.X(), a.Y() + b.Y(), a.Z() + b.Z()); }

inline
TMBVector3 operator - (const TMBVector3 &a, const TMBVector3 &b) {
   /// Subtraction of 3-vectors.
   return TVector3(a.X() - b.X(), a.Y() - b.Y(), a.Z() - b.Z()); }

inline
Double_t operator * (const TMBVector3 &a, const TMBVector3 &b) {
   /// Dot product of two vectors
   return a.Dot(b); }

inline
TMBVector3 operator * (const TMBVector3 &p, Double_t a) {
   /// Scalar product of 3-vectors.
   return TVector3(a*p.X(), a*p.Y(), a*p.Z()); }

inline
TMBVector3 operator * (Double_t a, const TMBVector3 &p) {
   /// Scaling of 3-vectors with a real number
   return TVector3(a*p.X(), a*p.Y(), a*p.Z()); }

inline
TMBVector3 operator * (const TMatrix &m, const TMBVector3 &v) {
   /// Transform with matrix
  return TVector3( m(0,0)*v.X()+m(0,1)*v.Y()+m(0,2)*v.Z(),
                   m(1,0)*v.X()+m(1,1)*v.Y()+m(1,2)*v.Z(),
                   m(2,0)*v.X()+m(2,1)*v.Y()+m(2,2)*v.Z());
}

/*************************************************************************
 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

// Author: Pasha Murat, Peter Malzacher   12/02/99

#endif // INCLUDE_TMBVECTOR3
Revision:	1.1
Committed:	Tue Nov 11 23:01:21 2008 UTC (16 years, 5 months ago) by kukartse
Content type:	text/plain
Branch:	MAIN
CVS Tags:	V00-03-01, ZMorph_BASE_20100408, gak040610_morphing, V00-02-02, gak011410, gak010310, ejterm2010_25nov2009, V00-02-01, V00-02-00, gak112409, CMSSW_22X_branch_base, segala101609, V00-01-15, V00-01-14, V00-01-13, V00-01-12, V00-01-11, V00-01-10, gak031009, gak030509, gak022309, gak021209, gak040209, gak012809, V00-01-09, V00-01-08, V00-01-07, V00-01-06, V00-01-05, V00-01-04, V00-00-07, V00-00-06, V00-00-05, V00-00-04, V00-01-03, V00-00-02, V00-00-01, HEAD
Branch point for:	ZMorph-V00-03-01, CMSSW_22X_branch
Log Message:	initial creation of a package for LJMET multivariate analysis
#	User	Rev	Content
1	kukartse	1.1	#ifndef INCLUDE_TMBVECTOR3
2			#define INCLUDE_TMBVECTOR3
3
4			#include "TMath.h"
5			#include "TError.h"
6			#include "TVector2.h"
7			#include "TMatrix.h"
8			#include "TRotation.h"
9
10			/*! \brief
11			The Physics Vector package
12
13			The Physics Vector package consists of five classes:
14			- TVector2
15			- TVector3
16			- TRotation
17			- TLorentzVector
18			- TLorentzRotation
19			It is a combination of CLHEPs Vector package written by
20			Leif Lonnblad, Andreas Nilsson and Evgueni Tcherniaev
21			and a ROOT package written by Pasha Murat.
22			for CLHEP see: http://wwwinfo.cern.ch/asd/lhc++/clhep/
23
24			\ingroup reco
25
26			<H2>
27			TVector3</H2>
28			<TT>TVector3</TT> is a general three vector class, which can be used for
29			the description of different vectors in 3D.
30			<H3>
31			Declaration / Access to the components</H3>
32			<TT>TVector3</TT> has been implemented as a vector of three <TT>Double_t</TT>
33			variables, representing the cartesian coordinates. By default all components
34			are initialized to zero:
35
36			<P><TT>  TVector3 v1;        //
37			v1 = (0,0,0)</TT>
38			<BR><TT>  TVector3 v2(1);     // v2 = (1,0,0)</TT>
39			<BR><TT>  TVector3 v3(1,2,3); // v3 = (1,2,3)</TT>
40			<BR><TT>  TVector3 v4(v2);    // v4 = v2</TT>
41
42			<P>It is also possible (but not recommended) to initialize a <TT>TVector3</TT>
43			with a <TT>Double_t</TT> or <TT>Float_t</TT> C array.
44
45			<P>You can get the basic components either by name or by index using <TT>operator()</TT>:
46
47			<P><TT>  xx = v1.X();    or    xx =
48			v1(0);</TT>
49			<BR><TT>  yy = v1.Y();
50			yy = v1(1);</TT>
51			<BR><TT>  zz = v1.Z();
52			zz = v1(2);</TT>
53
54			<P>The memberfunctions <TT>SetX()</TT>, <TT>SetY()</TT>, <TT>SetZ()</TT>
55			and<TT> SetXYZ()</TT> allow to set the components:
56
57			<P><TT>  v1.SetX(1.); v1.SetY(2.); v1.SetZ(3.);</TT>
58			<BR><TT>  v1.SetXYZ(1.,2.,3.);</TT>
59			<BR>
60			<H3>
61			Noncartesian coordinates</H3>
62			To get information on the <TT>TVector3</TT> in spherical (rho,phi,theta)
63			or cylindrical (z,r,theta) coordinates, the
64			<BR>the member functions <TT>Mag3()</TT> (=magnitude=rho in spherical coordinates),
65			<TT>Mag32()</TT>, <TT>Theta()</TT>, <TT>CosTheta()</TT>, <TT>Phi()</TT>,
66			<TT>Perp()</TT> (the transverse component=r in cylindrical coordinates),
67			<TT>Perp2()</TT> can be used:
68
69			<P><TT>  Double_t m  = v.Mag3();    // get magnitude
70			(=rho=Sqrt(xx+yy+z*z)))</TT>
71			<BR><TT>  Double_t m2 = v.Mag32();   // get magnitude squared</TT>
72			<BR><TT>  Double_t t  = v.Theta();  // get polar angle</TT>
73			<BR><TT>  Double_t ct = v.CosTheta();// get cos of theta</TT>
74			<BR><TT>  Double_t p  = v.Phi();    // get azimuth angle</TT>
75			<BR><TT>  Double_t pp = v.Perp();   // get transverse component</TT>
76			<BR><TT>  Double_t pp2= v.Perp2();  // get transvers component
77			squared</TT>
78
79			<P>It is also possible to get the transverse component with respect to
80			another vector:
81
82			<P><TT>  Double_t ppv1 = v.Perp(v1);</TT>
83			<BR><TT>  Double_t pp2v1 = v.Perp2(v1);</TT>
84
85			<P>The pseudorapiditiy ( eta=-ln (tan (phi/2)) ) can be get by <TT>Eta()</TT>
86			or <TT>PseudoRapidity()</TT>:
87			<BR>
88			<BR><TT>  Double_t eta = v.PseudoRapidity();</TT>
89
90			<P>There are set functions to change one of the noncartesian coordinates:
91
92			<P><TT>  v.SetTheta(.5); // keeping rho and phi</TT>
93			<BR><TT>  v.SetPhi(.8);   // keeping rho and theta</TT>
94			<BR><TT>  v.SetMag3(10.);  // keeping theta and phi</TT>
95			<BR><TT>  v.SetPerp(3.);  // keeping z and phi</TT>
96			<BR>
97			<H3>
98			Arithmetic / Comparison</H3>
99			The <TT>TVector3</TT> class provides the operators to add, subtract, scale and compare
100			vectors:
101
102			<P><TT>  v3  = -v1;</TT>
103			<BR><TT>  v1  = v2+v3;</TT>
104			<BR><TT>  v1 += v3;</TT>
105			<BR><TT>  v1  = v1 - v3</TT>
106			<BR><TT>  v1 -= v3;</TT>
107			<BR><TT>  v1 *= 10;</TT>
108			<BR><TT>  v1  = 5*v2;</TT>
109
110			<P><TT>  if(v1==v2) {...}</TT>
111			<BR><TT>  if(v1!=v2) {...}</TT>
112			<BR>
113			<H3>
114			Related Vectors</H3>
115			<TT>  v2 = v1.Unit();       // get unit
116			vector parallel to v1</TT>
117			<BR><TT>  v2 = v1.Orthogonal(); // get vector orthogonal to v1</TT>
118			<H3>
119			Scalar and vector products</H3>
120			<TT>  s = v1.Dot(v2);   // scalar product</TT>
121			<BR><TT>  s = v1 * v2;      // scalar product</TT>
122			<BR><TT>  v = v1.Cross(v2); // vector product</TT>
123			<H3>
124			Angle between two vectors</H3>
125			<TT>  Double_t a = v1.Angle(v2);</TT>
126			<H3>
127			Rotations</H3>
128
129			<H5>
130			Rotation around axes</H5>
131			<TT>  v.RotateX(.5);</TT>
132			<BR><TT>  v.RotateY(TMath::Pi());</TT>
133			<BR><TT>  v.RotateZ(angle);</TT>
134			<H5>
135			Rotation around a vector</H5>
136			<TT>  v1.Rotate(TMath::Pi()/4, v2); // rotation around v2</TT>
137			<H5>
138			Rotation by TRotation</H5>
139			<TT>TVector3</TT> objects can be rotated by objects of the <TT>TRotation</TT>
140			class using the <TT>Transform()</TT> member functions,
141			<BR>the <TT>operator =</TT> or the <TT>operator </TT> of the TRotation
142			class:
143
144			<P><TT>  TRotation m;</TT>
145			<BR><TT>  ...</TT>
146			<BR><TT>  v1.transform(m);</TT>
147			<BR><TT>  v1 = m*v1;</TT>
148			<BR><TT>  v1 = m; // Attention v1 = mv1</TT>
149			<H5>
150			Transformation from rotated frame</H5>
151			<TT>  TVector3 direction = v.Unit()</TT>
152			<BR><TT>  v1.RotateUz(direction); // direction must be TVector3 of
153			unit length</TT>
154
155			<P>transforms v1 from the rotated frame (z' parallel to direction, x' in
156			the theta plane and y' in the xy plane as well as perpendicular to the
157			theta plane) to the (x,y,z) frame.'
158
159			See copyright statements at end of file.
160
161			*/
162
163			class TMBVector3 : public TObject {
164			public:
165
166			TMBVector3(Double_t x = 0.0, Double_t y = 0.0, Double_t z = 0.0):
167			fX(x), fY(y), fZ(z) { /// The constructor.
168			}
169
170			TMBVector3(const Double_t *a): fX(a[0]), fY(a[1]), fZ(a[2])
171			{ /// Constructor from an array
172			}
173
174			TMBVector3(const Float_t *a): fX(a[0]), fY(a[1]), fZ(a[2])
175			{ /// Constructor from an array
176			}
177
178			TMBVector3(const TMBVector3 &p): TObject(p),
179			fX(p.fX), fY(p.fY), fZ(p.fZ)
180			{ /// Copy constructor.
181			}
182
183			TMBVector3(const TVector3& v): fX(v.X()), fY(v.Y()), fZ(v.Z())
184			{
185			/// conversion constructor from TVector3
186			}
187
188			virtual ~TMBVector3() {}
189			// Destructor
190
191			Double_t operator () (int i) const {
192			/// Get components by index.
193			switch(i) {
194			case 0: return fX;
195			case 1: return fY;
196			case 2: return fZ;
197			default: Error("operator()(i)", "bad index (%d) returning &fX",i);
198			}
199			return 0.;
200			}
201
202			inline Double_t operator [] (int i) const { /// Get components by index.
203			return operator()(i); }
204
205
206			Double_t & operator () (int i) {
207			/// Get components by index.
208
209			switch(i) {
210			case 0: return fX;
211			case 1: return fY;
212			case 2: return fZ;
213			default: Error("operator()(i)", "bad index (%d) returning &fX",i);
214			}
215			return fX;
216			}
217			inline Double_t & operator [] (int i) { /// Set components by index.
218			return operator()(i); }
219
220			Double_t x() const { /// get X component
221			return fX; }
222			Double_t X() const { /// get X component
223			return x(); }
224			Double_t Px() const { /// get X component
225			return x(); }
226
227			Double_t y() const { /// get Y component
228			return fY; }
229			Double_t Y() const { /// get Y component
230			return y(); }
231			Double_t Py() const { /// get X component
232			return y(); }
233
234			Double_t z() const { /// get Y component
235			return fZ; }
236			Double_t Z() const { /// get Z component
237			return z(); }
238			Double_t Pz() const { /// get X component
239			return z(); }
240
241			void SetX(Double_t a) { /// set X component, see TVector3::SetX()
242			fX=a; }
243			void SetPx(Double_t a) { /// set X component, see TVector3::SetX()
244			SetX(a); }
245			void SetY(Double_t a) { /// set X component, see TVector3::SetY()
246			fY=a; }
247			void SetPy(Double_t a) { /// set X component, see TVector3::SetX()
248			SetY(a); }
249			void SetZ(Double_t a) { /// set X component, see TVector3::SetZ()
250			fZ=a; }
251			void SetPz(Double_t a) { /// set X component, see TVector3::SetX()
252			SetZ(a); }
253			inline void SetXYZ(Double_t x, Double_t y, Double_t z) {
254			/// set all members using Cartesian coordinates
255			fX=x; fY=y; fZ=z; }
256
257			void SetPtEtaPhi(Double_t pt, Double_t eta, Double_t phi) {
258			/// set all members using \|momentum\|, eta, and phi
259			Double_t apt = TMath::Abs(pt);
260			Double_t d=TMath::Tan(2.0*TMath::ATan(TMath::Exp(-eta)));
261			if (d==0.)
262			Warning("SetPtEtaPhi", "tan(2*atan(exp(-eta)))==0.! Setting z=0.");
263			SetXYZ(aptTMath::Cos(phi), aptTMath::Sin(phi), d==0.?0.:apt/d );
264			}
265			void SetPtThetaPhi(Double_t pt, Double_t theta, Double_t phi) {
266			/// set all members using \|momentum\|, theta, and phi
267			Double_t apt = TMath::Abs(pt);
268			fX = apt * TMath::Cos(phi);
269			fY = apt * TMath::Sin(phi);
270			Double_t tanTheta = TMath::Tan(theta);
271			if (tanTheta==0.)
272			Warning("SetPtThetaPhi", "tan(theta)==0.! Setting z=0.");
273			fZ = tanTheta ? apt / tanTheta : 0;
274			}
275
276			inline void GetXYZ(Double_t *carray) const {
277			/// Getters into an arry (no checking of array bounds!)
278			carray[0]=fX; carray[1]=fY; carray[2]=fZ; }
279			inline void GetXYZ(Float_t *carray) const {
280			/// Getters into an arry (no checking of array bounds!)
281			carray[0]=fX; carray[1]=fY; carray[2]=fZ; }
282
283			Double_t Eta() const { /// get eta (aka pseudo-rapidity)
284			return PseudoRapidity(); }
285			Double_t Theta() const { /// get theta
286			return fX==.0 && fY==.0 && fZ==.0 ? .0 : TMath::ATan2(Perp(),fZ); }
287			Double_t CosTheta() const { /// get cos(Theta)
288			Double_t ptot = Mag3();
289			return ptot == 0.0 ? 1.0 : fZ/ptot;
290			}
291			Double_t Phi() const { /// get phi from 0..2Pi
292			return (fX==.0 && fY==.0 ? .0 : TVector2::Phi_0_2pi(TMath::ATan2(fY,fX))); }
293
294			Double_t Rho() const { /// get magnitude of momentum ("radius" in spherical coords)
295			return Mag3(); }
296
297			virtual Double_t Mag32() const {
298			/// The magnitude squared (rho^2 in spherical coordinate system).
299			return fXfX + fYfY + fZ*fZ; }
300			virtual Double_t Mag3() const {
301			/// The magnitude (rho in spherical coordinate system).
302			return TMath::Sqrt(Mag32()); }
303			Double_t P() const { /// get \|P\|
304			return Mag3(); }
305
306			Double_t DeltaPhi(const TMBVector3 &v) const {
307			/// Get difference in phi, between -pi and pi
308			return TVector2::Phi_mpi_pi(Phi()-v.Phi()); }
309
310			Double_t DeltaR(const TMBVector3 &v) const {
311			/// Get "distance" dR of v to this, dR=sqrt(dPhidPhi+dEta*dEta)
312			Double_t deta = Eta()-v.Eta();
313			Double_t dphi = TVector2::Phi_mpi_pi(Phi()-v.Phi());
314			return TMath::Sqrt( detadeta+dphidphi );
315			}
316			Double_t DrEtaPhi(const TMBVector3 &lv) const {
317			/// Get "distance" dR of lv to this, dR=sqrt(dPhidPhi+dEta*dEta)
318			return DeltaR(lv); }
319
320			// TVector2 EtaPhiVector() { /// return TVector2(eta,phi)
321			// return TVector2 (Eta(),Phi()); }
322
323			Double_t Angle(const TVector3 & q) const { /// Angle wrt. another vector.
324			Double_t ptot2 = Mag32()*q.Mag2();
325			if(ptot2 <= 0.) return .0;
326			Double_t arg = Dot(q)/TMath::Sqrt(ptot2);
327			if(arg > 1.0) arg = 1.0;
328			if(arg < -1.0) arg = -1.0;
329			return TMath::ACos(arg);
330			}
331
332			void SetPhi(Double_t ph) { /// set phi keeping mag and theta constant
333			Double_t xy = Perp();
334			SetX(xy*TMath::Cos(ph));
335			SetY(xy*TMath::Sin(ph));
336			}
337
338			inline void SetTheta(Double_t th) { /// Set theta keeping mag and phi constant
339			Double_t ma = Mag3();
340			Double_t ph = Phi();
341			SetX(maTMath::Sin(th)TMath::Cos(ph));
342			SetY(maTMath::Sin(th)TMath::Sin(ph));
343			SetZ(ma*TMath::Cos(th));
344			}
345
346			inline void SetMag3(Double_t ma) { /// Set magnitude keeping theta and phi constant
347			Double_t factor = Mag3();
348			if (factor == 0) {
349			Warning("SetMag3","zero vector can't be stretched");
350			}else{
351			factor = ma/factor;
352			SetX(fX*factor);
353			SetY(fY*factor);
354			SetZ(fZ*factor);
355			}
356			}
357			inline Double_t Perp2() const {
358			/// The transverse component squared (R^2 in cylindrical coordinate system)
359			return fXfX + fYfY; }
360			inline Double_t Pt() const {
361			/// transverse component; projection of 3 vector onto XY plane (R in cylindrical coords)
362			return Perp(); }
363			inline Double_t Perp() const {
364			/// The transverse component (R in cylindrical coordinate system)
365			return TMath::Sqrt(Perp2()); }
366
367			inline void SetPerp(Double_t r) {
368			/// Set the transverse component keeping phi and z constant.
369			Double_t p = Perp();
370			if (p != 0.0) {
371			fX *= r/p;
372			fY *= r/p;
373			}
374			}
375
376			inline Double_t Perp2(const TMBVector3 &p) const {
377			/// The transverse component w.r.t. given axis squared.
378			Double_t tot = p.Mag32();
379			Double_t ss = Dot(p);
380			return tot > 0.0 ? Mag32()-ss*ss/tot : Mag32();
381			}
382
383			inline Double_t Pt(const TMBVector3 &p) const {
384			/// The transverse component w.r.t. given axis.
385			return Perp(p); }
386			inline Double_t Perp(const TMBVector3 &p) const {
387			/// The transverse component w.r.t. given axis.
388			return TMath::Sqrt(Perp2(p)); }
389
390			inline void SetMag3ThetaPhi(Double_t mag, Double_t theta, Double_t phi) {
391			/// Set all components as magnitude, theta, and phi
392			Double_t amag = TMath::Abs(mag);
393			fX = amag * TMath::Sin(theta) * TMath::Cos(phi);
394			fY = amag * TMath::Sin(theta) * TMath::Sin(phi);
395			fZ = amag * TMath::Cos(theta);
396			}
397
398			inline TMBVector3 & operator = (const TMBVector3 &p) {
399			/// Assignment.
400			fX = p.fX; fY = p.fY; fZ = p.fZ;
401			return *this;
402			}
403
404			inline Bool_t operator == (const TMBVector3 &v) const {
405			/// Equality operator. Equal if all components X,Y,Z are equal.
406			return (v.fX==fX && v.fY==fY && v.fZ==fZ); }
407			inline Bool_t operator != (const TMBVector3 &v) const {
408			/// Inequality operator. Not equal if any of X,Y,Z are not equal.
409			return (v.fX!=fX \|\| v.fY!=fY \|\| v.fZ!=fZ); }
410
411			/// Equivalence operator. True if all components are
412			/// equivalent within machine precision.
413			Bool_t is_equal (const TMBVector3 &v) const;
414
415			inline TMBVector3 & operator += (const TMBVector3 &p) { /// Addition.
416			fX += p.fX; fY += p.fY; fZ += p.fZ;
417			return *this;
418			}
419
420			inline TMBVector3 & operator -= (const TMBVector3 &p) { /// Subtraction.
421			fX -= p.fX; fY -= p.fY; fZ -= p.fZ;
422			return *this;
423			}
424
425			inline TMBVector3 operator - () const { /// Unary minus.
426			return TMBVector3(-fX, -fY, -fZ); }
427
428			inline TMBVector3 & operator *= (Double_t a) {
429			/// Scaling with real numbers.
430			fX = a; fY = a; fZ *= a;
431			return *this;
432			}
433
434			inline TMBVector3 Unit() const { /// Unit vector parallel to this.
435			Double_t tot = Mag32();
436			TMBVector3 p(fX,fY,fZ);
437			return tot>.0 ? p *= (1.0/TMath::Sqrt(tot)) : p;
438			}
439
440			inline TMBVector3 Orthogonal() const { /// Vector orthogonal to this.
441			Double_t x = fX < 0.0 ? -fX : fX;
442			Double_t y = fY < 0.0 ? -fY : fY;
443			Double_t z = fZ < 0.0 ? -fZ : fZ;
444			if (x < y)
445			return x<z ? TMBVector3(0,fZ,-fY) : TMBVector3(fY,-fX,0);
446			return y<z ? TMBVector3(-fZ,0,fX) : TMBVector3(fY,-fX,0);
447			}
448
449			Double_t Dot(const TMBVector3 &p) const { /// Scalar product.
450			return fXp.fX + fYp.fY + fZ*p.fZ; }
451
452			inline TMBVector3 Cross(const TMBVector3 &p) const { /// Cross product
453			return TMBVector3(fYp.fZ-p.fYfZ, fZp.fX-p.fZfX, fXp.fY-p.fXfY);
454			}
455
456			Double_t PseudoRapidity() const {
457			/// Returns the pseudo-rapidity, i.e. -ln(tan(theta/2))
458			double cosTheta = CosTheta();
459			if (cosTheta*cosTheta < 1)
460			return (-0.5* TMath::Log( (1.0-cosTheta)/(1.0+cosTheta) ));
461			Warning("PseudoRapidity","transvers momentum = 0! return +/- 10e10");
462			if (fZ > 0) return (10e10);
463			else return (-10e10);
464			}
465			void RotateX(Double_t angle) {
466			/// Rotates around the x-axis.
467			Double_t s = TMath::Sin(angle);
468			Double_t c = TMath::Cos(angle);
469			Double_t yy = fY;
470			fY = cyy - sfZ;
471			fZ = syy + cfZ;
472			}
473
474			void RotateY(Double_t angle) {
475			/// Rotates around the y-axis.
476			Double_t s = TMath::Sin(angle);
477			Double_t c = TMath::Cos(angle);
478			Double_t zz = fZ;
479			fZ = czz - sfX;
480			fX = szz + cfX;
481			}
482
483			void RotateZ(Double_t angle) {
484			/// Rotates around the z-axis.
485			Double_t s = TMath::Sin(angle);
486			Double_t c = TMath::Cos(angle);
487			Double_t xx = fX;
488			fX = cxx - sfY;
489			fY = sxx + cfY;
490			}
491
492			void RotateUz(const TMBVector3& NewUzVector) {
493			/// Rotates reference frame from Uz to newUz (must be unit vector!)
494			Double_t u1 = NewUzVector.fX;
495			Double_t u2 = NewUzVector.fY;
496			Double_t u3 = NewUzVector.fZ;
497			Double_t up = u1u1 + u2u2;
498
499			if (up>0.) {
500			up = TMath::Sqrt(up);
501			Double_t px = fX, py = fY, pz = fZ;
502			fX = (u1u3px - u2py + u1up*pz)/up;
503			fY = (u2u3px + u1py + u2up*pz)/up;
504			fZ = (u3u3px - px + u3uppz)/up;
505			}
506			else if (u3 < 0.) { fX = -fX; fZ = -fZ; } // phi=0 teta=pi
507			else {};
508			}
509
510			void Rotate(Double_t angle, const TMBVector3 &axis) {
511			/// Rotates around the axis specified by another vector
512			TRotation trans;
513			trans.Rotate(angle, (TVector3)axis);
514			operator*=(trans);
515			}
516
517			TMBVector3 & operator *= (const TRotation &m) {
518			/// Transform with a rotation matrix
519			return this = m (*this);
520			}
521			TMBVector3 & Transform(const TRotation &m) {
522			/// Transform with a rotation matrix
523			return this = m (*this);
524			}
525
526			TVector2 XYvector() const {
527			/// TVector2 containing x and y components
528			return TVector2(fX,fY); }
529
530			operator TVector3() const {
531			/// cast to a TVector3
532			return TVector3(X(),Y(),Z()); }
533
534			void Print(Option_t* option="") const {
535			/// print components to stdout
536			Printf("%s %s (x,y,z)=(%f,%f,%f) (rho,theta,phi)=(%f,%f,%f)",
537			GetName(),GetTitle(),X(),Y(),Z(),
538			Mag3(),Theta()TMath::RadToDeg(),Phi()TMath::RadToDeg());
539			}
540
541			// Helper to test if two FP numbers are equivalent within machine precision.
542			static bool is_equal (double x1, double x2);
543
544			private:
545			Double32_t fX; ///< X component (left-right)
546			Double32_t fY; ///< Y component (up-down)
547			Double32_t fZ; ///< Z component (along the beam)
548
549			//ClassDef(TMBVector3,1) // A 3D physics vector
550
551			};
552
553			inline
554			TMBVector3 operator + (const TMBVector3 &a, const TMBVector3 &b) {
555			/// Addition of 3-vectors.
556			return TVector3(a.X() + b.X(), a.Y() + b.Y(), a.Z() + b.Z()); }
557
558			inline
559			TMBVector3 operator - (const TMBVector3 &a, const TMBVector3 &b) {
560			/// Subtraction of 3-vectors.
561			return TVector3(a.X() - b.X(), a.Y() - b.Y(), a.Z() - b.Z()); }
562
563			inline
564			Double_t operator * (const TMBVector3 &a, const TMBVector3 &b) {
565			/// Dot product of two vectors
566			return a.Dot(b); }
567
568			inline
569			TMBVector3 operator * (const TMBVector3 &p, Double_t a) {
570			/// Scalar product of 3-vectors.
571			return TVector3(ap.X(), ap.Y(), a*p.Z()); }
572
573			inline
574			TMBVector3 operator * (Double_t a, const TMBVector3 &p) {
575			/// Scaling of 3-vectors with a real number
576			return TVector3(ap.X(), ap.Y(), a*p.Z()); }
577
578			inline
579			TMBVector3 operator * (const TMatrix &m, const TMBVector3 &v) {
580			/// Transform with matrix
581			return TVector3( m(0,0)v.X()+m(0,1)v.Y()+m(0,2)*v.Z(),
582			m(1,0)v.X()+m(1,1)v.Y()+m(1,2)*v.Z(),
583			m(2,0)v.X()+m(2,1)v.Y()+m(2,2)*v.Z());
584			}
585
586			/*************************************************************************
587			* Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
588			* All rights reserved. *
589			* *
590			* For the licensing terms see $ROOTSYS/LICENSE. *
591			* For the list of contributors see $ROOTSYS/README/CREDITS. *
592			*************************************************************************/
593
594			// Author: Pasha Murat, Peter Malzacher 12/02/99
595
596			#endif // INCLUDE_TMBVECTOR3