// sv-5.cc
// 		A ZOOM PhysicsVectors tutorial example
//
// 		Tutorial example 5 for the SpaceVector and UnitVector classes:
//
// 		Projections and Angles

// * Projection along a direction
// * Perpendicular parts
// * The angle between two vectors
// * Decomposition into polar and azimuthal parts

// Key lines are marked with 					// <-----------

	// Get the SpaceVector and UnitVector classes

#include "ZMutility/ZMenvironment.h"
#include "ZMutility/FixedTypes.h"
#include "PhysicsVectors/SpaceVector.h"
#include "PhysicsVectors/UnitVector.h"
USING( std::cout )
USING( std::endl )

	// This example will use iostream output, so get that --
	// in a portable way.

#include "ZMutility/iostream"
using std::cout;
using std::endl;

const Float8 PI =  3.1415926535897931;

int main() {

	// (See tutorial examples sv-1 and sv-3 for how to construct
	//  and access UnitVectors and SpaceVectors.)

	// generate a few sample vectors
  SpaceVector v1 (  3,  2,  5 );
  SpaceVector v2 (  4, -2,  3 );

  UnitVector u1 ( 0, RADIANS, PI/2, RADIANS );
  UnitVector u2 ( 45, DEGREES, PI, RADIANS );


// ***** Projection along a direction

	// Here is how to obtain a projection along the Z axis
	//----------------------------------------------------

  SpaceVector vz = v1.project();				// <-----------

  cout << "Projection of " << v1 << " along Z is " << vz << endl;

	// Here is how to obtain a projection along an arbitrary axis
	//-----------------------------------------------------------

  SpaceVector vu1 = v1.project(u1);				// <-----------

  cout << "Projection of " << v1 << "\n along " << u1 << " is " << vz << endl;

  SpaceVector v3 ( 3, -2, 1 );
  SpaceVector vv3 = v1.project(v3);				// <-----------
	// This projects along the direction defined by the normalized v3
	// Note that this need not be the same as v3.project(v1)


// ***** Perpendicular parts

	// Here is how to obtain scalar perp() or its square
	//--------------------------------------------------

  Float8 p = v2.perp();						// <----------
  Float8 p2 = v2.perp2();					// <----------

	// perp2() can be done much quicker than perp()

  p  = v2.perp(u2);						// <----------
  p2 = v2.perp2(u2);						// <----------

	// Here is how to obtain vector perpendicular part
	//------------------------------------------------

  SpaceVector  vtrans = v1.perpPart();				// <----------

  vz = v1.project();
  if ( v1.isNear(vtrans+vz) ) {
    cout << v1 << " == transverse + longitudinal parts\n" <<
	vtrans << " + " << vz << endl;
  }

  SpaceVector  vtransU1 = v1.perpPart(u1);			// <----------

  SpaceVector vlongU1 = v1.project(u1);
  if ( v1.isNear(vtransU1+vlongU1) ) {
    cout << v1 << " == transverse + longitudinal parts\n" <<
	vtrans << " + " << vz << endl;
  }

// ***** The angle between two vectors

	// Here is how to find theta or cos theta
	//---------------------------------------

  Float8 thet = v1.angle(v2);					// <----------

  cout << "Angle between " << v1 << " and " << v2 << " is " << thet << endl;

  Float8 cthet = v1.cosTheta(v2);				// <----------
	// cosTheta() is quicker

  cout << "The easier-to-calculate cosine of that is " << cthet << endl;

        // Here is how to use a UnitVector as the reference direction
        //-----------------------------------------------------------

  Float8 angl = v1.theta(u1);                                   // <----------
        // theta() is a synonym for angle()

  cout << "Angle between " << v1 << " and " << v2 << " is " << angl << endl;

	// Here is how to find eta, the pseudorapidity along a direction
	//--------------------------------------------------------------

  Float8 pr = v1.eta(v2);					// <----------
  	// - log tan (v1.angle(v2)/2)

  cout << "Psuedorapidity of " << v1 << "\n in direction specified by "
	<< UnitVector(v2) << " is " << pr << endl;

// ***** Decomposition into polar and azimuthal parts

	// Here is how to decompose an angle relative to the Z axis
	//---------------------------------------------------------

  Float8 deltheta = v1.polarAngle(v2);				// <----------
	// difference in thetas
  Float8 delphi   = v1.azimAngle(v2);				// <----------
	// difference in phi

  cout << "Angle between " << v1 << " and " << v2 << "\ndecomposes into polar "
	<< deltheta << " and azimuthal " << delphi << endl;

	// Here is how to decompose an angle relative to an arbitrary axis
	//----------------------------------------------------------------

  deltheta = v1.polarAngle(v2,u1);				// <----------
	// difference in angles against u1
  delphi   = v1.azimAngle(v2,u1);				// <----------
	// angle between projections into plane normal to u1

  cout << "Angle between " << v1 << " and " << v2 << " about " << u1
       << "\ndecomposes into polar " << deltheta << " and azimuthal "
       << delphi << endl;

  return 0;

} /* end of main */