// lv-1.cc
// 		A ZOOM PhysicsVectors tutorial example
//
// 		Tutorial example 1 for the LorentzVector class:
//
// 		Construction and component access
//
// * Constructing LorentzVectors and output to ostreams
// * Changing individual components
// * Working with the SpaceVector part of a LorentzVector
// * X_HAT4, Y_HAT4, Z_HAT4, and T_HAT4

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

	// Get the LorentzVector class -- this also pulls in SpaceVector

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

        // These classes themselves pull in FixedTypes, but just to
        // remind ourselves that we are using Float8, not double, for
        // 64-bit floats...

#include "ZMutility/FixedTypes.h"
ZM_USING_NAMESPACE( zmfxt )

	// 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() {

// ***** Constructing LorentzVectors and Output to ostreams

	// Here is how to construct a 4-vector from x, y, z, t
	//----------------------------------------------------

  Float8 px = 2;
  Float8 py = 3;
  Float8 pz = 25;
  Float8 e  = 30;
  LorentzVector w1 ( px, py, pz, Tcomponent(e) );	// <------------------

  cout << "w1 is " << w1 << endl;

	// Here is how to construct a 4-vector from a SpaceVector and t
	//-------------------------------------------------------------

  SpaceVector v1 ( px, py, pz );
  LorentzVector w2 ( v1, e );			 	// <-----------------

  	// The way to construct a LorentzVector whose space part
	// is expressed in polar coordinates is to first make the
	// space part:

  SpaceVector v2 ( 10, 45, DEGREES, 30, DEGREES );
  LorentzVector w3 ( v2, e );
  cout << "w3 is " << w3 << endl;

  	// You can do this in one step:

  LorentzVector w4 (SpaceVector(5, 2.5, ETA, 60, RADIANS), 30); // <---------
  cout << "w4 is " << w4 << endl;

// * Constructing an array of LorentzVectors

	// Here is how to construct an array of  4-vectors
	//-------------------------------------------------

  LorentzVector www[100];				// <------------------
	// These are all initialized to (0,0,0,0)


// ***** Changing individual components

	// Here is how to change one Cartesian component, keeping others fixed
	//--------------------------------------------------------------------

  w1.setX( 3 );						// <------------------
  w1.setT( 13 );					// <------------------

  w1.setY( w1.x() - w1.t() );				// <------------------

  w1.setZ( w1.z() * .2 );				// <------------------

  cout << "w1 = " << w1 << endl;

	// Here is how to change just the Space part keeping t fixed
	//----------------------------------------------------------

  w1.v() = SpaceVector ( 2, 90, DEGREES, PI/8, RADIANS ); // <---------------

  cout << "w1 = " << w1 << endl;

	// Here is how to change all four components at once
	//--------------------------------------------------

  w3.set ( 3, 4, 5, Tcomponent(10) );			// <-----------------

  w3.set ( SpaceVector(10, 90, DEGREES, 0, DEGREES ), 20); // <--------------


// ***** Working with the SpaceVector part of a LorentzVector

	// Here is how to determine polar or cylindrical coordinates
	//----------------------------------------------------------

  w1.v() = SpaceVector ( 4, 30, DEGREES, 15, DEGREES );

  Float8 phi   = w1.getV().phi();			// <-----------------
  Float8 theta = w1.getV().theta();			// <-----------------
  Float8 r     = w1.getV().r();				// <-----------------
  Float8 rho   = w1.getV().rho();			// <-----------------

  cout << "phi = " << phi << " theta = " << theta
	<< " rho = " << rho << " r = " << r << endl;

	// The following syntax will also be supported:
	//  phi = w1.v().phi();

	// However, w1.v() is not itself a true SpaceVector object:
 	// To take the angle of that against v1, for instance,
	// you could not say w1.v().angle(v1) -- instead, you do:

  Float8 proj = w1.getV().angle(v1);			// <-----------------

// ***** X_HAT4, Y_HAT4, Z_HAT4, and T_HAT4

	// The package supplies 4-vectors of unit length along each axis:

  cout << "X_HAT4 is " << X_HAT4 << endl;		// <-----------------

	// in contrast to:

  cout << "X_HAT  is " << X_HAT  << endl;

  w1 = 5 * X_HAT4 - 3 * Y_HAT4 + Z_HAT4 + 10 * T_HAT4;	// <-----------------

  cout << "w1 is " << w1 << endl;

  return 0;

} /* end of main */