Vec.hh

Go to the documentation of this file.

#ifndef VEC_HH_INCLUDED
#define VEC_HH_INCLUDED

#include <cmath>
#include <iostream>
#include <sstream>
#include "Rtypes.h"

/**
 * The Vec class is a straightforward 3-d vector, which also works in pyroot.
 */
struct Vec
{
  double x,y,z;

  /**
   * Constructor.
   *
   * \param  x_            x position
   * \param  y_            y position
   * \param  z_            z position
   */
  Vec(double x_, double y_, double z_) : x(x_), y(y_), z(z_) {}
  
  /**
   * Default constructor.
   */
  Vec():x(0),y(0),z(0) {}
  
  /**
   * Get dot product.
   *
   * \param  v             vector
   * \return               dot product
   */
  double dot(const Vec& v) const { return v.x*x + v.y*y+ v.z*z;}
  
  /**
   * Get cross product.
   *
   * \param  r             vector
   * \return               cross product
   */
  Vec cross(const Vec r)   const { return Vec ( y*r.z-z*r.y, z*r.x-x*r.z, x*r.y-y*r.x);}

  /**
   * Add vector.
   *
   * \param  v             vector
   * \return               this vector
   */
  Vec& operator+=(const Vec& v)  { x+=v.x; y+=v.y; z+=v.z; return *this;}

  /**
   * Subtract vector.
   *
   * \param  v             vector
   * \return               this vector
   */
  Vec& operator-=(const Vec& v)  { x-=v.x; y-=v.y; z-=v.z; return *this;}

  /**
   * Multiply vector.
   *
   * \param  d             factor
   * \return               this vector
   */
  Vec& operator*=(double d)      { x*=d; y*=d; z*=d; return *this;}

  /**
   * Divide vector.
   *
   * \param  d             factor
   * \return               this vector
   */
  Vec& operator/=(double d)      { return operator*=( 1.0 / d ); }

  /**
   * Negate vector.
   *
   * \return               vector
   */
  Vec operator-() const          { return Vec(-x,-y,-z); }
  
  /**
   * Check equality with given vector.
   *
   * \param  v             vector
   * \return               true if (x,y,z) positions of two vectors are equal; else false
   */
  bool operator==( const Vec& v ) const { return x==v.x && y==v.y && z==v.z ; }

  /**
   * Check in-equality with given vector.
   *
   * \param  v             vector
   * \return               true if one of (x,y,z) positions of two vectors are not equal; else false
   */
  bool operator!=( const Vec& v ) const { return x!=v.x || y!=v.y || z!=v.z ; }

  /**
   * Set vector.
   *
   * \param  xx            x position
   * \param  yy            y position
   * \param  zz            z position
   * \return               this vector
   */
  Vec& set(double xx, double yy, double zz) { x=xx; y=yy; z=zz; return *this;}

  /**
   * Set vector according given zenith and azimuth angles.
   *
   * \param  theta         zenith  angle [rad]
   * \param  phi           azimuth angle [rad]
   * \return               this vector
   */
  Vec& set_angles(double theta, double phi)
  {
    x = sin ( theta ) * cos( phi );
    y = sin ( theta ) * sin( phi );
    z = cos ( theta );
    return *this;
  }

  /**
   * Get azimuth angle.
   *
   * \return               angle [rad]
   */
  double phi()   const { return atan2( y,x ); }

  /**
   * Get zenith angle.
   *
   * \return               angle [rad]
   */
  double theta() const { return acos(z); }

  /**
   * Get length.
   *
   * \return               length
   */
  double len()   const { double l = dot(*this); return (l > 0)? sqrt(l) : 0; }

  /**
   * Get length of (x,y) component.
   *
   * \return               length
   */
  double lenxy() const { const double r2 = x*x + y*y; return (r2>0) ? sqrt(r2) :0; }

  /**
   * Normalise this vector.
   *
   * \return               this vector
   */
  Vec& normalize() { return operator/=( len() ); }

  /**
   * Print vector.
   *
   * \param  out           output stream
   */
  void print( std::ostream& out = std::cout ) const
  {
        out << "Vec:" << x << " " << y << " " << z;
  }

  /**
   * Get string representation of this vector
   *
   * \return               string
   */
  const char* __repr__() const 
  {
    static std::string buffer;

    std::ostringstream s; 
  
    print(s);

    buffer = s.str();

    return buffer.c_str();
  }

  /**
   * Add vector.
   *
   * \param  v             vector
   * \return               vector
   */
  Vec __add__(const Vec& v) const { Vec r=*this; return r+=v; }

  /**
   * Subtract vector.
   *
   * \param  v             vector
   * \return               vector
   */
  Vec __sub__(const Vec& v) const { Vec r=*this; return r-=v; }

  /**
   * Multiply vector.
   *
   * \param  d             factor
   * \return               vector
   */
  Vec __mul__(double d )    const { Vec r=*this; return r*=d; }

  /**
   * Multiply vector.
   *
   * \param  d             factor
   * \return               vector
   */
  Vec __rmul__(double d )   const { return __mul__(d);        } 

  /**
   * Divide vector.
   *
   * \param  d             factor
   * \return               vector
   */
  Vec __div__(double d )    const { Vec r=*this; return r/=d; }

  /**
   * Rotate around z-axis with given angle.
   *
   * \param  ang           angle [rad]
   * \return               this vector
   */
  Vec& rotate_z(double ang)
  {
    const Vec o = *this; 
    x = o.x *cos(ang) - o.y * sin(ang);
    y = o.x *sin(ang) + o.y * cos(ang);
    z = o.z;
    return *this;
  }

  /**
   * Rotate around x-axis with given angle.
   *
   * \param  ang           angle [rad]
   * \return               this vector
   */
  Vec& rotate_x(double ang)
  {
    const Vec o = *this;
    x = o.x;
    y = o.y *cos(ang) + o.z * -sin(ang);
    z = o.y *sin(ang) + o.z * cos(ang);
    return *this;
  }

  /**
   * Rotate around y-axis with given angle.
   *
   * \param  ang           angle [rad]
   * \return               this vector
   */
  Vec& rotate_y(double ang)
  {
    const Vec o = *this;
    x = o.x *cos(ang) + o.z * sin(ang);
    y = o.y;
    z = -o.x *sin(ang) + o.z * cos(ang);
    return *this;
  }

  ClassDefNV(Vec,3)
};

/**
 * Write vector to output stream.
 *
 * \param  out           output stream
 * \param  v             vector
 * \return               output stream
 */
inline std::ostream& operator<<( std::ostream& out , const Vec& v )
{
  out << v.x << " " << v.y << " " << v.z << " ";
  return out;
}

/**
 * Read vector from input stream.
 *
 * \param  in            input stream
 * \param  v             vector
 * \return               input stream
 */
inline std::istream& operator>>(std::istream& in, Vec& v) 
{ 
  in >> v.x >> v.y >> v.z ; return in;
}

/**
 * Get cosine of space angle between two vectors.
 *
 * \param  a             first  vector
 * \param  b             second vector
 * \return               cosine
 */
inline double cos_angle_between( const Vec& a, const Vec& b)
{
  const double n = a.len() * b.len();
  return a.dot(b) / n;
}

/**
 * Get space angle between two vectors.
 *
 * \param  a             first  vector
 * \param  b             second vector
 * \return               angle [rad]
 */
inline double angle_between( const Vec& a, const Vec& b )
{
  double c = cos_angle_between( a, b );
  if ( c < -1 ) return M_PI;
  if ( c > 1  ) return 0;
  return acos( c );
}

/**
 * Add two vectors.
 *
 * \param  a             first  vector
 * \param  b             second vector
 * \return               vector
 */
inline Vec operator+(const Vec& a, const Vec& b) { Vec r(a); return r+=b;}

/**
 * Subtract two vectors.
 *
 * \param  a             first  vector
 * \param  b             second vector
 * \return               vector
 */
inline Vec operator-(const Vec& a, const Vec& b) { Vec r(a); return r-=b;}

/**
 * Multiply vector.
 *
 * \param  a             factor
 * \param  v             vector
 * \return               vector
 */
inline Vec operator*(double a, const Vec& v)     { return Vec(a*v.x,a*v.y,a*v.z);}

/**
 * Multiply vector.
 *
 * \param  v             vector
 * \param  a             factor
 * \return               vector
 */
inline Vec operator*(const Vec& v, double a)     { return Vec(a*v.x,a*v.y,a*v.z);}

/**
 * Divide vector.
 *
 * \param  v             vector
 * \param  a             factor
 * \return               vector
 */
inline Vec operator/(const Vec& v, double a)     { return Vec(v.x/a,v.y/a,v.z/a);}

#endif