// sv-4.cc
// 		A ZOOM PhysicsVectors tutorial example
//
// 		Tutorial example 4 for the SpaceVector class:
//
// 		Comparing Vectors
//
// * Exact Comparison
// * Relative difference
// * DeltaR
// * Parallel and Orthogonal

// 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"
ZM_USING_NAMESPACE( zmpv )
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;

int main() {

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

	// The methods illustrated here may all be applied to UnitVector
	// as well.  (See example sv-3 for how to construct UnitVectors.)

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


// ***** Exact Comparison

	// Here is how to apply exact comparisons to vectors
	//--------------------------------------------------

  v2 = v1;
  v3 = v2 + v4;
  if ( v2 == v1 ) { 					// <------------------
    cout << v2 << " == " << v1 << endl;
  }
  if ( v4 != v1 ) { 					// <------------------
    cout << v4 << " != " << v1 << endl;
  }

  	// Other comparisons > < >= <= are also provided so  that
	// sort templates can be used on lists of SpaceVectors.


// ***** Relative Difference

	// Here is how to determine if two vectors are "nearly equal"
	//-----------------------------------------------------------

  v1.set (  1, 2, 3  );
  v3.set ( -3, 3, 1  );
  Float8 epsilon = 1.0E-17;
  v2 = v1 + epsilon * v3;	// Some calculations have generated two
  				// vectors which are effectively equal:
  if ( v2 != v1 ) { 		// These vectors are not exactly equal
    cout << v2 << " != " << v1 << endl;
  }

  if ( v2.isNear (v1) ) { 				// <------------------
    cout << "   but " << v2 << ".isNear" << v1 << endl;
  }

  	// isNear does a relative (not absolute) comparison:
	// Re-scaling does not affect the result.

  SpaceVector v1BIG = v1 * 1.0E20;
  SpaceVector v2BIG = v2 * 1.0E20;
  if ( v2BIG.isNear (v1BIG) ) { 			// <------------------
    cout << "   but " << v2BIG << ".isNear" << v1BIG <<endl;
  }

	// isNear compares the magnitude of the vector difference,
	// normalized to the sqrt of the dot product, to some small
	// tolerance (100 times the machine epsilon).  You can also
	// just get that ratio:

  Float8 x = v1.howNear(v2);				// <------------------
  cout << v1 << ".howNear" << v2 << " = " << x << endl;
  cout << v1BIG << ".howNear" << v2BIG << " = "<< v1BIG.howNear(v2BIG) << endl;

	// You can also get the square of the absolute difference of
	// two vectors.  This of course is not scale invariant.

  x = v1BIG.diff2(v2BIG);				// <------------------
  cout << v1BIG << ".diff2" << v2BIG << " = "<< x << endl;



	// Here is how to control the tolerance used to determine near equality
	//---------------------------------------------------------------------

  v2 = v1 + .0001 * v3;

  if ( v2.isNear(v1, .002) ) {				// <------------------
    cout << "Within a tolerance of .002, "
	   << v2 << ".isNear" << v1 << endl;
  }

  	// Just to prove the supplied tolerance was used:
  Float8 t = SpaceVector::getTolerance();		// <------------------
  if ( ! v2.isNear(v1) ) {
    cout << "Within the default tolerance of " << t
	   << ",\n" << v2 << " is not near " << v1 << endl;
  }

	// The default tolerance used is 2.2E-16 but you may specify a
	// different tolerance as above.  Or you may change the default:

  SpaceVector::setTolerance(.002);			// <------------------
  if ( v2.isNear(v1) ) {
    cout << "Within the new default tolerance of "
	 << SpaceVector::getTolerance()
	 << ",\n" << v2 << ".isNear" << v1 << endl;
  }

  SpaceVector::setTolerance(t);
	// I like to leave things the way I found them...


// ***** DeltaR

	// Here is how to find deltaR = sqrt ((delta_eta)**2 + (delta_phi)**2)
	//--------------------------------------------------------------------

  v1.set ( .05, .07, 125.0 );
  v2.set ( .04, .10, 125.1 );

  Float8 d = v1.deltaR(v2);				// <------------------

  cout << "DeltaR between " << v1 << " and " << v2 << " is " << d << endl;

// ***** Parallel and Orthogonal

	// Here is how to determine if two vectors are nearly parallel
	//------------------------------------------------------------

  v1.set ( 4, 135, DEGREES, -15, DEGREES );
  epsilon = 1.0E-17;
  v2 = v1 + epsilon * UnitVector ( 1, 1, 0 );
	// These vectors are not exactly parallel

  if ( v2.isParallel(v1) ) { 	    			// <------------------
    cout << v2 << ".isParallel" << v1 << endl;
  }

 	// The tolerance used is the default tolerance, but you can override it:

  v2 = v1 + .01 * UnitVector ( 1, 3, 0 );

  if ( v2.isParallel(v1, .2) ) {			// <------------------
    cout << "Within a tolerance of .2, "
	<< v2 << ".isParallel" << v1 << endl;
  }

	// And you can find the measure used to decide parallel or not:
  	// This is the ratio of the cross product to the dot product.

  Float8 y = v1.howParallel(v2);			// <------------------
  cout << v1 << ".howParallel" << v2 << " = " << y << endl;


	// Here is how to determine if two vectors are nearly orthogonal
	//--------------------------------------------------------------

  v1.set ( 4.3,  6, -8  );
  v2.set ( -4,   5, 5.2 );
  v3 = v1.cross(v2);       // This is orthogonal - to machine accuracy - to v1

  if ( v1.isOrthogonal(v3) ) { 	    			// <------------------
    cout << v1 << ".isOrthogonal" << v3 << endl;
  }

 	// The tolerance used is the default tolerance, but you can override it:

  v1 += .01 * SpaceVector(-3, 4, 2) ;

  if ( v3.isOrthogonal(v1, .1) ) {			// <------------------
    cout << "Within a tolerance of .1, "
	<< v3 << ".isOrthogonal" << v1 << endl;
  }

	// And you can find the measure used to decide parallel or not:
  	// This is the ratio of the dot product to the cross product.
	// All these relative nearness measures are pegged to a maximum of 1.

  y = v1.howOrthogonal(v3);				// <------------------
  cout << v1 << ".howOrthogonal" << v3 << " = " << y << endl;

  return 0;

} /* end of main */