// -*- C++ -*-
//
// This file is a side header, dependant on CLHEP.
//
// It declares, defines, and implements the UnitVector class present in ZOOM
// PhysicsVectors.  It does so with no non-inline implementation code; thus
// no code goes into the CLHEP library, and if this header is not included, 
// the existence of UnitVector is quite invisible.
//
// Most methods merely call the corresponding Hep3Vector method.  
// The mechanism for doing this in non-trivial cases is a C-style cast 
//	WHICH DEPENDS CRITICALLY ON THE FACT THAT THE REPRESENTATIONS OF
//	Hep#Vector and UnitVector MATCH!
//
//	IF THE PRIVATE DATA OF Hep3Vector CHANGES, THE PRIVATE DATA HERE MUST!
//
// One could use the safer operator Hep3Vector, as in 
// Hepdouble getTheta const {return (this->operator Hep3Vector()).gettheta();}
// but this implies creation of a temporary and in some cases that would risk
// a significant performance hit.  (Optimization can trivially eliminate that
// hit but can we trust optimizers these days?)
//

// .SS Authors
// Mark Fischler
//

#ifndef UNITVECTOR_H
#define UNITVECTOR_H

#ifdef GNUPRAGMA
#pragma interface
#endif

#include "CLHEP/config/CLHEP.h"
#include "PhysicsVectors/SpaceVector.h"

// NOTE:  UnitVector will be typedefed to HepUnit3Vector inside the 
//	  optional zmpv namespace.  Thus if and when HepUnit3Vector is 
//	  accepted into CLHEP, the relationship will be analogous to that
//	  between SpaceVector and Hep3Vector.


class HepUnit3Vector {

protected:  // The IDENTICAL variables as Hep3Vector, in the same order
  HepDouble dx;
  HepDouble dy;
  HepDouble dz;
  
public:

  // **** Conversion to and construction from Hep3Vector ****

  operator SpaceVector() const {return SpaceVector( dx, dy, dz );}
  operator Hep3Vector() const {return Hep3Vector( dx, dy, dz );}
  // Thus when used as a const ARGUMENT, 
  // a HepUnit3Vector may replace a SpaceVector

  // ****Constructors and setting entire vector ****

  HepUnit3Vector() : dx(0.0), dy(0.0), dz(1.0) {}

  void set   (HepDouble x, HepDouble y, HepDouble z) 
    {((Hep3Vector*)this)->set(x, y, z); ((Hep3Vector*)this)->setR(1.0);}
  HepUnit3Vector (HepDouble x, HepDouble y, HepDouble z) : dx(x), dy(y), dz(z) 
    {((Hep3Vector*)this)->setR(1.0);}
  // Normalize these cartesian coordinates.

  void set   (const Hep3Vector & v) 
    {((Hep3Vector*)this)->set(v.x(), v.y(), v.z()); 
     ((Hep3Vector*)this)->setR(1.0);}
  HepUnit3Vector (const Hep3Vector & v) : dx(v.x()), dy(v.y()), dz(v.z()) 
    {((Hep3Vector*)this)->setR(1.0);}
  // Normalize the cartesian coordinates defined by this vector.

  HepUnit3Vector & rectify()
    {HepDouble radius= sqrt( dx*dx+dy*dy+dz*dz);
     dx=dx/radius; dy=dy/radius; dz=dz/radius; return *this;}

  // No cylindrical coordinates

  HepUnit3Vector (const HepUnit3Vector &u) : dx(u.dx), dy(u.dy), dz(u.dz) {}
  // The copy constructor.

  ~HepUnit3Vector() {}
  // The destructor.  Not virtual - inheritance from this class is dangerous.

  // **** Component Access ****

  HepDouble operator () (int i) const 
		{return ((Hep3Vector*)this)->operator()(i);}
  HepDouble operator [] (int i) const 
		{return ((Hep3Vector*)this)->operator[](i);}
  // Get components by index 

  // No Setting components by index.

  HepDouble getX() const {return dx;}
  HepDouble getY() const {return dy;}
  HepDouble getZ() const {return dz;}
  // The components in cartesian coordinate system. 

  HepDouble x() const {return dx;}
  HepDouble y() const {return dy;}
  HepDouble z() const {return dz;}
  // The components in cartesian coordinate system.  Same as getX() etc.

  // No Setting the components in cartesian coordinate system.
  // No Setting all three components in cartesian coordinate system.

  HepDouble phi()    const {return ((Hep3Vector*)this)->phi();}
  HepDouble getPhi() const {return phi();}
  // The azimuth angle.

  void setPhi(HepDouble x) { ((Hep3Vector*)this)->setPhi(x);}
  // Set phi keeping theta constant 

  HepDouble theta()    const {return ((Hep3Vector*)this)->theta();}
  HepDouble getTheta() const {return theta();}
  HepDouble angle()    const {return theta();}
  // The polar angle.

  void setTheta(HepDouble x) { ((Hep3Vector*)this)->setTheta(x);} 
  // Set theta keeping phi constant 

  HepDouble cosTheta() const {return dz;}
  // Cosine of the polar angle.

  HepDouble cos2Theta() const {return dz*dz;} 
  // Cosine squared of the polar angle.

  HepDouble mag2() const {return 1.0;}
  // The magnitude squared (r^2 in spherical coordinate system).

  HepDouble mag()  const {return 1.0;}
  HepDouble getR() const {return 1.0;}
  HepDouble r()    const {return 1.0;}
  // The magnitude (r in spherical coordinate system).

  // No setMag() or setR()

  inline void setR ( HepDouble s );
  // setMag()

  HepDouble perp2() const { return dx*dx + dy*dy; }
  // The transverse component squared (rho^2 in cylindrical coordinate system).

  HepDouble perp() const { return sqrt(perp2()); }
  inline HepDouble rho    () const {return perp();} 
  inline HepDouble getRho () const {return perp();}
  // The transverse component (rho in cylindrical coordinate system).

  // No setPerp or setRho
  // No setCylTheta
  // Set theta while keeping transvers component and phi fixed 

  HepDouble perp2(const Hep3Vector & v) const 
		{return ((Hep3Vector*)this)->perp2(v);}
  // The transverse component w.r.t. given axis squared.

  HepDouble perp2(const HepUnit3Vector & u) const  
	{ HepDouble uv = dot(u);
	  HepDouble p2 = 1.0 - uv*uv;
	  return p2>0 ? p2:0; }

  HepDouble perp(const Hep3Vector & v) const {return sqrt(perp2(v));}
  // The transverse component w.r.t. given axis.

  HepDouble perp(const HepUnit3Vector & u) const {return sqrt(perp2(u));}

  // **** Assignment **** 

  HepUnit3Vector & operator = (const HepUnit3Vector & u) 
		{dx=u.dx; dy=u.dy; dz=u.dz; return *this;}

  // NO operator += 
  // NO operator -= 
  // NO operator *= 
  // NO operator /= 

  HepUnit3Vector operator - () const {return HepUnit3Vector(-dx, -dy, -dz);}
  // Unary minus.

  // **** Comparisions ****

  HepInt compare (const Hep3Vector & v) const
		{return ((Hep3Vector*)this)->compare(v);}

  HepBoolean operator == (const Hep3Vector & v) const
		{return ((Hep3Vector*)this)->operator==(v);}
  HepBoolean operator != (const Hep3Vector & v) const
		{return ((Hep3Vector*)this)->operator!=(v);}

  HepBoolean operator > (const Hep3Vector & v) const
		{return ((Hep3Vector*)this)->operator>(v);}
  HepBoolean operator < (const Hep3Vector & v) const
		{return ((Hep3Vector*)this)->operator<(v);}
  HepBoolean operator>= (const Hep3Vector & v) const
		{return ((Hep3Vector*)this)->operator>=(v);}
  HepBoolean operator<= (const Hep3Vector & v) const
		{return ((Hep3Vector*)this)->operator<=(v);}
  // dictionary ordering according to z, then y, then x component

  HepBoolean isNear (const Hep3Vector & v, 
		     HepDouble epsilon=Hep3Vector::getTolerance()) const
		{return ((Hep3Vector*)this)->isNear(v,epsilon);}

  HepBoolean isNear (const HepUnit3Vector & u, 
		     HepDouble epsilon=Hep3Vector::getTolerance()) const {
	return (*this-u).mag2() <= epsilon*epsilon;
  }

  HepDouble howNear(const Hep3Vector & v ) const {return (*this - v).mag();}

  HepDouble deltaR(const Hep3Vector & v) const 
		{return ((Hep3Vector*)this)->deltaR(v);}
  
  inline HepDouble diff2 (const Hep3Vector & v) const
		{return ((Hep3Vector*)this)->diff2(v);}
  // |v1-v2|**2

  HepBoolean isParallel (const Hep3Vector & v, 
		     HepDouble epsilon=Hep3Vector::getTolerance()) const
		{return ((Hep3Vector*)this)->isParallel(v,epsilon);}
  // Are the vectors parallel, within the given tolerance?

  HepBoolean isOrthogonal (const Hep3Vector & v, 
		     HepDouble epsilon=Hep3Vector::getTolerance()) const
		{return ((Hep3Vector*)this)->isOrthogonal(v,epsilon);}
  // Are the vectors orthogonal, within the given tolerance?

  HepDouble howParallel   (const Hep3Vector & v) const
		{return ((Hep3Vector*)this)->howParallel(v);}
  // | v1.cross(v2) / v1.dot(v2) |, to a maximum of 1.

  HepDouble howOrthogonal (const Hep3Vector & v) const
		{return ((Hep3Vector*)this)->howOrthogonal(v);}
  // | v1.dot(v2) / v1.cross(v2) |, to a maximum of 1.

  // No setTolerance or getTolerance; uses the Hep3Vector tolerance

  // **** Intrinsic properties ****

  Hep3Vector unit() const {return (Hep3Vector)(*this);}

  // No orthogonal()
  // No beta()
  // No gamma() 
  // No coLinearRapidity() 
  // No rapidity()

  HepDouble pseudoRapidity() const
		{return ((Hep3Vector*)this)->pseudoRapidity();}
  HepDouble eta     () const {return pseudoRapidity();}
  HepDouble getEta  () const {return pseudoRapidity();}

  void setEta  ( HepDouble p )
		{((Hep3Vector*)this)->setEta(p);}

  // NO setCylEta  

  // **** Combining with other vectors ****

  HepDouble dot(const Hep3Vector & v) const 
		{return ((Hep3Vector*)this)->dot(v);}

  Hep3Vector cross(const Hep3Vector & v) const
		{return ((Hep3Vector*)this)->cross(v);}

  HepDouble angle(const Hep3Vector & v) const
		{return ((Hep3Vector*)this)->angle(v);}
  HepDouble theta(const Hep3Vector & v) const	{return angle(v);}

  HepDouble cosTheta (const HepUnit3Vector & u) const {return dot(u);}
  HepDouble cosTheta (const Hep3Vector & v) const 
		{return ((Hep3Vector*)this)->cosTheta(v);}

  HepDouble cos2Theta(const HepUnit3Vector & u) const
		{HepDouble c = dot(u); return c*c;}
  HepDouble cos2Theta(const Hep3Vector & v) const 
		{return ((Hep3Vector*)this)->cos2Theta(v);}
  // cos^2 of the angle between two vectors

  Hep3Vector project () const {return Hep3Vector(0,0,dz);}
  Hep3Vector project (const HepUnit3Vector & u) const
		{return dot(u)*u;}
  Hep3Vector project (const Hep3Vector & v) const
		{return ((Hep3Vector*)this)->project(v);}
  // projection of a vector along a direction.  

  Hep3Vector perpPart() const {return Hep3Vector(dx,dy,0);}
  Hep3Vector perpPart (const HepUnit3Vector & u) const
		{return (*this)-dot(u)*u;}
  Hep3Vector perpPart (const Hep3Vector & v) const
		{return ((Hep3Vector*)this)->perpPart(v);}
  // vector minus its projection along a direction.

  // no rapidity (v)

  HepDouble eta(const Hep3Vector & v) const 
		{return ((Hep3Vector*)this)->eta(v);}
  // - ln tan of the angle beween the vector and the ref direction.

  // **** Angle decomposition ****
  // Note:  All of these could be done faster, especially using a UnitVector
  //	    as the reference.  However, these are not extensively used and
  //	    thus we use the simplest possible coding, based on Hep3Vector.

  HepDouble polarAngle (const Hep3Vector & v2) const
		{return ((Hep3Vector*)this)->polarAngle(v2);}
  // The reference direction is Z: the polarAngle is abs(u.theta()-v2.theta()).

  HepDouble azimAngle  (const Hep3Vector & v2) const
		{return ((Hep3Vector*)this)->azimAngle(v2);}
  // The reference direction is Z: the azimAngle is u.phi()-v2.phi()

  HepDouble polarAngle (const Hep3Vector & v2, 
					const Hep3Vector & ref) const
		{return ((Hep3Vector*)this)->polarAngle(v2,ref);}
  // For arbitrary reference direction, 
  // 	polarAngle is abs(v.angle(ref) - v2.angle(ref)).

  HepDouble azimAngle  (const Hep3Vector & v2, 
					const Hep3Vector & ref) const
		{return ((Hep3Vector*)this)->azimAngle(v2,ref);}
  // To compute azimangle, project u and v2 into the plane normal to
  // the reference direction.  Then in that plane take the angle going
  // clockwise around the direction from projection of u to that of v2.

  // **** Rotating ****

  Hep3Vector & rotateX(HepDouble delta)
		{return ((Hep3Vector*)this)->rotateX(delta);}

  Hep3Vector & rotateY(HepDouble delta)
		{return ((Hep3Vector*)this)->rotateY(delta);}

  Hep3Vector & rotateZ(HepDouble delta)
		{return ((Hep3Vector*)this)->rotateZ(delta);}

  // No rotateUz(const Hep3Vector&);

  Hep3Vector & rotate (const Hep3Vector & axis, HepDouble delta)
		{return ((Hep3Vector*)this)->rotate(axis,delta);}
  Hep3Vector & rotate (HepDouble delta, const Hep3Vector & v) 
		{return rotate(v,delta);}

  // no operator *= (const HepRotation &);
  // no transform(const HepRotation &);

  Hep3Vector & rotate  (const HepAxisAngle & ax)
		{return ((Hep3Vector*)this)->rotate(ax);}

  Hep3Vector & rotate (const HepEulerAngles & e)
		{return ((Hep3Vector*)this)->rotate(e);}

  Hep3Vector & rotate (HepDouble phi,
                        HepDouble theta,
                        HepDouble psi)
		{return ((Hep3Vector*)this)->rotate(phi,theta,psi);}
  // Rotate via Euler Angles. Our Euler Angles conventions are 
  // those of Goldstein Classical Mechanics page 107.

  // **** Global methods friend declatations ****

  inline friend HepUnit3Vector rotationXOf (const HepUnit3Vector & u, HepDouble delta);
  inline friend HepUnit3Vector rotationYOf (const HepUnit3Vector & u, HepDouble delta);
  inline friend HepUnit3Vector rotationZOf (const HepUnit3Vector & u, HepDouble delta);
  inline friend HepUnit3Vector rotationOf  (const HepUnit3Vector & u, 
			const Hep3Vector & axis, HepDouble delta);
  inline friend HepUnit3Vector rotationOf 
			(const HepUnit3Vector & u, const HepAxisAngle & ax);
  inline friend HepUnit3Vector rotationOf (const HepUnit3Vector & u, 
			HepDouble phi, HepDouble theta, HepDouble psi);
  inline friend HepUnit3Vector rotationOf 
			(const HepUnit3Vector & u, const HepEulerAngles & e);

};


// **** Global Methods ****

HepUnit3Vector rotationXOf (const HepUnit3Vector & u, HepDouble delta) {	
	Hep3Vector v( rotationXOf((Hep3Vector)u,delta) );
	HepUnit3Vector w; w.dx = v.x(); w.dy = v.y(); w.dz=v.z(); return w;}

HepUnit3Vector rotationYOf (const HepUnit3Vector & u, HepDouble delta) {
	Hep3Vector v( rotationYOf((Hep3Vector)u,delta) );
	HepUnit3Vector w; w.dx = v.x(); w.dy = v.y(); w.dz=v.z(); return w;}

HepUnit3Vector rotationZOf (const HepUnit3Vector & u, HepDouble delta) {
	Hep3Vector v( rotationZOf((Hep3Vector)u,delta) );
	HepUnit3Vector w; w.dx = v.x(); w.dy = v.y(); w.dz=v.z(); return w;}

HepUnit3Vector rotationOf (const HepUnit3Vector & u, 
				const Hep3Vector & axis, HepDouble delta) {
	Hep3Vector v( rotationOf((Hep3Vector)u,axis,delta) );
	HepUnit3Vector w; w.dx = v.x(); w.dy = v.y(); w.dz=v.z(); return w;}

HepUnit3Vector rotationOf (const HepUnit3Vector & u, const HepAxisAngle & ax) {
	Hep3Vector v( rotationOf((Hep3Vector)u,ax) );
	HepUnit3Vector w; w.dx = v.x(); w.dy = v.y(); w.dz=v.z(); return w;}

HepUnit3Vector rotationOf (const HepUnit3Vector & u, 
			HepDouble phi, HepDouble theta, HepDouble psi) {
	Hep3Vector v( rotationOf((Hep3Vector)u,phi,theta,psi) );
	HepUnit3Vector w; w.dx = v.x(); w.dy = v.y(); w.dz=v.z(); return w;}

HepUnit3Vector rotationOf (const HepUnit3Vector & u, const HepEulerAngles & e) {
	Hep3Vector v( rotationOf((Hep3Vector)u,e) );
	HepUnit3Vector w; w.dx = v.x(); w.dy = v.y(); w.dz=v.z(); return w;}

// operator << (HepStd::ostream & os, const Hep3Vector &u) 
//	not needed since there is a conversion to Hep3Vector

inline HepStd::istream & operator >> (HepStd::istream &is, HepUnit3Vector &u) {
  Hep3Vector v;
  is >> v;
  u.set (v.x(), v.y(), v.z());
  return is;
}
// Input from a stream, normalizing the three numbers.

ZM_BEGIN_NAMESPACE( zmpv )      /*  namespace zmpv  {  */

typedef double Scalar;

class UnitVector : public HepUnit3Vector { 

public:

  // Constructors parallel to those for HepUnit3Vector

  UnitVector() { dx = 0.0; dy = 0.0; dz = 1.0; }

  UnitVector (HepDouble x, HepDouble y, HepDouble z) 
    {((Hep3Vector*)this)->set(x,y,z); ((Hep3Vector*)this)->setR(1.0);}
  void set   (HepDouble x, HepDouble y, HepDouble z) 
    {((Hep3Vector*)this)->set(x, y, z); ((Hep3Vector*)this)->setR(1.0);}
  // Normalize these cartesian coordinates.

  UnitVector (const Hep3Vector & v) 
    {((Hep3Vector*)this)->operator=(v); ((Hep3Vector*)this)->setR(1.0);}
  // Normalize the cartesian coordinates defined by this vector.

  UnitVector (const HepUnit3Vector &u) : HepUnit3Vector(u) {}
  UnitVector (const UnitVector &u)     : HepUnit3Vector(u) {}
  // Copy constructors.

  ~UnitVector() {}
  // The destructor.  Not virtual - inheritance from this class is dangerous.

  // Spherical coordinate constructors:

  void set   (HepDouble theta, ZMpvRADIANS_t thetaUnits,
              HepDouble phi,   ZMpvRADIANS_t phiUnits) 
    {((SpaceVector*)this)->set(1.0, theta, thetaUnits, phi, phiUnits);}
  void set   (HepDouble theta, ZMpvRADIANS_t thetaUnits,
              HepDouble phi,   ZMpvDEGREES_t phiUnits) 
    {((SpaceVector*)this)->set(1.0, theta, thetaUnits, phi, phiUnits);}
  void set   (HepDouble theta, ZMpvDEGREES_t thetaUnits,
              HepDouble phi,   ZMpvRADIANS_t phiUnits) 
    {((SpaceVector*)this)->set(1.0, theta, thetaUnits, phi, phiUnits);}
  void set   (HepDouble theta, ZMpvDEGREES_t thetaUnits,
              HepDouble phi,   ZMpvDEGREES_t phiUnits) 
    {((SpaceVector*)this)->set(1.0, theta, thetaUnits, phi, phiUnits);}
  void set   (HepDouble theta, ZMpvETA_t     thetaUnits,
              HepDouble phi,   ZMpvRADIANS_t phiUnits) 
    {((SpaceVector*)this)->set(1.0, theta, thetaUnits, phi, phiUnits);}
  void set   (HepDouble theta, ZMpvETA_t     thetaUnits,
              HepDouble phi,   ZMpvDEGREES_t phiUnits) 
    {((SpaceVector*)this)->set(1.0, theta, thetaUnits, phi, phiUnits);}

  UnitVector (HepDouble theta, ZMpvRADIANS_t thetaUnits,
                  HepDouble phi,   ZMpvRADIANS_t phiUnits) 
	    {set(theta, thetaUnits, phi, phiUnits);}
  UnitVector (HepDouble theta, ZMpvRADIANS_t thetaUnits,
                  HepDouble phi,   ZMpvDEGREES_t phiUnits) 
	    {set(theta, thetaUnits, phi, phiUnits);}
  UnitVector (HepDouble theta, ZMpvDEGREES_t thetaUnits,
                  HepDouble phi,   ZMpvRADIANS_t phiUnits) 
	    {set(theta, thetaUnits, phi, phiUnits);}
  UnitVector (HepDouble theta, ZMpvDEGREES_t thetaUnits,
                  HepDouble phi,   ZMpvDEGREES_t phiUnits) 
	    {set(theta, thetaUnits, phi, phiUnits);}
  UnitVector (HepDouble theta, ZMpvETA_t     thetaUnits,
                  HepDouble phi,   ZMpvRADIANS_t phiUnits) 
	    {set(theta, thetaUnits, phi, phiUnits);}
  UnitVector (HepDouble theta, ZMpvETA_t     thetaUnits,
                  HepDouble phi,   ZMpvDEGREES_t phiUnits) 
	    {set(theta, thetaUnits, phi, phiUnits);}

}; // UnitVector	

// **** Axis UnitVectors ****

static const UnitVector X_HAT (1,0,0);
static const UnitVector Y_HAT (0,1,0);
static const UnitVector Z_HAT (0,0,1);

ZM_END_NAMESPACE( zmpv )        /*  }  // namespace zmpv  */

#endif /* UNITVECTOR_H */