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
#	Content
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