TMBVector3.hpp

Go to the documentation of this file.

#ifndef INCLUDE_TMBVECTOR3
#define INCLUDE_TMBVECTOR3

#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) {  
   }

   TMBVector3(const Double_t *a): fX(a[0]), fY(a[1]), fZ(a[2]) 
   { 
   }

   TMBVector3(const Float_t *a): fX(a[0]), fY(a[1]), fZ(a[2]) 
   { 
   }

   TMBVector3(const TMBVector3 &p): TObject(p),
      fX(p.fX), fY(p.fY), fZ(p.fZ)
   {  
   }

   TMBVector3(const TVector3& v): fX(v.X()), fY(v.Y()), fZ(v.Z())
   {
   }

  virtual ~TMBVector3() {}
  // Destructor

  Double_t operator () (int i) const {
     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 { 
     return operator()(i); }


  Double_t & operator () (int i) {

     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) { 
     return operator()(i); }

   Double_t x() const { 
      return fX; }
   Double_t X() const { 
      return x(); }
   Double_t Px() const { 
      return x(); }

   Double_t y() const { 
      return fY; }
   Double_t Y() const { 
      return y(); }
   Double_t Py() const { 
      return y(); }

   Double_t z() const { 
      return fZ; }
   Double_t Z() const { 
      return z(); }
   Double_t Pz() const { 
      return z(); }

   void SetX(Double_t a) { 
      fX=a; }
   void SetPx(Double_t a) { 
      SetX(a); }
   void SetY(Double_t a) { 
      fY=a; }
   void SetPy(Double_t a) { 
      SetY(a); }
   void SetZ(Double_t a) { 
      fZ=a; }
   void SetPz(Double_t a) { 
      SetZ(a); }
   inline void SetXYZ(Double_t x, Double_t y, Double_t z) {
       fX=x; fY=y; fZ=z; }
    
   void SetPtEtaPhi(Double_t pt, Double_t eta, Double_t 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) {
      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 {
      carray[0]=fX; carray[1]=fY; carray[2]=fZ; }
   inline void GetXYZ(Float_t *carray) const {
      carray[0]=fX; carray[1]=fY; carray[2]=fZ; }

   Double_t Eta() const { 
      return PseudoRapidity(); }
   Double_t Theta() const { 
      return fX==.0 && fY==.0 && fZ==.0 ? .0 : TMath::ATan2(Perp(),fZ); }
   Double_t CosTheta() const { 
      Double_t ptot = Mag3();
      return ptot == 0.0 ? 1.0 : fZ/ptot;
   }
   Double_t Phi() const { 
      return (fX==.0 && fY==.0 ? .0 : TVector2::Phi_0_2pi(TMath::ATan2(fY,fX))); }

   Double_t Rho() const { 
      return Mag3(); }

   virtual Double_t Mag32() const {
      return fX*fX + fY*fY + fZ*fZ; }
   virtual Double_t Mag3() const {
      return TMath::Sqrt(Mag32()); }
   Double_t P()  const { 
      return Mag3(); }

   Double_t DeltaPhi(const TMBVector3 &v) const { 
      return TVector2::Phi_mpi_pi(Phi()-v.Phi()); }

   Double_t DeltaR(const TMBVector3 &v) const {
      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 {
      return DeltaR(lv); }

    //   TVector2 EtaPhiVector() { /// return TVector2(eta,phi)
    //      return TVector2 (Eta(),Phi()); }

   Double_t Angle(const TVector3 & q) const { 
      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) { 
      Double_t xy   = Perp();
      SetX(xy*TMath::Cos(ph));
      SetY(xy*TMath::Sin(ph));
   }

  inline void SetTheta(Double_t th) { 
     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) { 
     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 { 
     return fX*fX + fY*fY; }
  inline Double_t Pt() const { 
     return Perp(); }
  inline Double_t Perp() const { 
     return TMath::Sqrt(Perp2()); }
     
  inline void SetPerp(Double_t r) { 
     Double_t p = Perp();
     if (p != 0.0) {
        fX *= r/p;
        fY *= r/p;
     }
  }

  inline Double_t Perp2(const TMBVector3 &p) const { 
     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 {
     return Perp(p); }
  inline Double_t Perp(const TMBVector3 &p) const {
     return TMath::Sqrt(Perp2(p)); }

  inline void SetMag3ThetaPhi(Double_t mag, Double_t theta, Double_t 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) {
     fX = p.fX; fY = p.fY; fZ = p.fZ;
     return *this;
  }

  inline Bool_t operator == (const TMBVector3 &v) const {
     return (v.fX==fX && v.fY==fY && v.fZ==fZ); }
  inline Bool_t operator != (const TMBVector3 &v) const {
     return (v.fX!=fX || v.fY!=fY || v.fZ!=fZ); }

  Bool_t is_equal (const TMBVector3 &v) const;

  inline TMBVector3 & operator += (const TMBVector3 &p) { 
     fX += p.fX; fY += p.fY; fZ += p.fZ;
     return *this;
  }

  inline TMBVector3 & operator -= (const TMBVector3 &p) { 
     fX -= p.fX; fY -= p.fY; fZ -= p.fZ;
     return *this;
  }

  inline TMBVector3 operator - () const { 
     return TMBVector3(-fX, -fY, -fZ); }

  inline TMBVector3 & operator *= (Double_t a) { 
     fX *= a; fY *= a; fZ *= a;
     return *this;
  }

  inline TMBVector3 Unit() const { 
     Double_t  tot = Mag32();
     TMBVector3 p(fX,fY,fZ);
     return tot>.0 ? p *= (1.0/TMath::Sqrt(tot)) : p;
  }

  inline TMBVector3 Orthogonal() const { 
     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 { 
     return fX*p.fX + fY*p.fY + fZ*p.fZ; }

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

  Double_t PseudoRapidity() const {
     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) {
     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) {
     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) {
     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) {
     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) {
     TRotation trans;
     trans.Rotate(angle, (TVector3)axis);
     operator*=(trans);
  }

  TMBVector3 & operator *= (const TRotation &m) {
     return *this = m * (*this);
  }
  TMBVector3 & Transform(const TRotation &m) {
     return *this = m * (*this);
  }

   TVector2 XYvector() const {
     return TVector2(fX,fY); }

   operator TVector3() const { 
      return TVector3(X(),Y(),Z()); }

   void Print(Option_t* option="") const {
      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; 
   Double32_t fY; 
   Double32_t fZ; 

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

};

inline
TMBVector3 operator + (const TMBVector3 &a, const TMBVector3 &b) {
  return TVector3(a.X() + b.X(), a.Y() + b.Y(), a.Z() + b.Z()); }

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

inline
Double_t operator * (const TMBVector3 &a, const TMBVector3 &b) {
   return a.Dot(b); }

inline
TMBVector3 operator * (const TMBVector3 &p, Double_t a) {
   return TVector3(a*p.X(), a*p.Y(), a*p.Z()); }

inline
TMBVector3 operator * (Double_t a, const TMBVector3 &p) {
   return TVector3(a*p.X(), a*p.Y(), a*p.Z()); }

inline
TMBVector3 operator * (const TMatrix &m, const TMBVector3 &v) {
  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

---