// lv-2.cc
// 		A ZOOM PhysicsVectors tutorial example
//
// 		Tutorial example 2 for the LorentzVector class:
//
// 		Properties and Arithmetic
// * Linear arithmetic
// * Dot product
// * Single-4-vector properties
// * Various norms and comparisons

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

	// Get the LorentzVector class

#include "ZMutility/ZMenvironment.h"
#include "ZMutility/FixedTypes.h"
#include "PhysicsVectors/LorentzVector.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 lv-1 for how to construct and access
	//  LorentzVectors.)

	// generate a few sample vectors

  const LorentzVector w1 (  3,  2,  5, Tcomponent (12) );
  const LorentzVector w2 (  4, -2,  3, Tcomponent (9) );
  LorentzVector w3 (  1,  2, -6, Tcomponent (14) );
  LorentzVector w4 ( -1,  3,  4, Tcomponent (8) );

  LorentzVector w5, w6, w7;

// ***** Linear arithmetic

	// Here is how to add, subtract, multiply by scalars, and so forth
	//----------------------------------------------------------------

  w5 = w1 + w2 - w3;					   // <--------------
  cout << w1 << " + " << w2
	<< "\n- " << w3 << " = " << w5 << endl;

  w6 = w1 * 2 + 3 * w2;					   // <--------------
  w7 = w7 / 2;						   // <--------------
  w7 += w5;						   // <--------------
  cout << "w7 = " << w7 << endl;

// ***** Dot product

	// Here is how to take a dot product w1.dot(w2)
	//---------------------------------------------

  Float8 d = w1.dot(w2);				   // <--------------
  cout << w1 << ".dot" << w2 << " = " << d << endl;

	// The dot product is metric-dependant:

  LorentzVector::setMetric (TimeNegative);
  d = w1.dot(w2);
  cout << "With (+ + + -) metric, " << w1 << ".dot"
	<< w2 << " = " << d << endl;
  LorentzVector::setMetric (TimePositive);

// ***** Single-4-vector properties

	// Here is how to find the invariant magnitude squared
	//----------------------------------------------------

  Float8 m2 = w1.mag2();				   // <--------------
  cout << w1 << ".mag2() = " << m2 << endl;

	// This is equivalent to w1.dot(w1) and is metric-dependant

	// Here is how to find the Euclidean Norm
	//---------------------------------------

  Float8 eNorm2 = w2.euclideanNorm2();			   // <--------------
  cout << w1 << ".euclideanNorm2() = " << eNorm2 << endl;

	// This is x**2 + y**2 + z**2 + t**2

  cout << "euclideanNorm2 + mag2 = " << eNorm2 + m2
	<< " = 2 * t**2 = " << 2*w1.t()*w1.t() << endl;
  	// But with TimeNegative metric,
	// eNorm2 - m2 = 2*t**2

  Float8 eNorm = w2.euclideanNorm();			   // <--------------

	// Here is how to find plus and minus parts
	//-----------------------------------------

  Float8 wplus  = w5.plus();				   // <--------------
  Float8 wminus = w5.minus();				   // <--------------

	// These are relative to the z axis:  t + z, t - z

  cout << w5 << ".plus() = " << wplus << endl;
  cout << w5 << ".minus() = " << wplus << endl;

	// Less useful for HEP are plus, minus relative to
	// an arbitrary axis:

  Float8 wXplus  = w5.plus(X_HAT);			   // <--------------
  Float8 wXminus = w5.minus(X_HAT);			   // <--------------

  cout << w5 << ".plus" << X_HAT << " = " << wXplus << endl;
  cout << w5 << ".minus" << X_HAT << " = " << wXminus << endl;

  UnitVector u ( 45, DEGREES, 0, DEGREES );
	// UnitVector and SpaceVector are pulled in by
	// LorentzVector.h

  Float8 wUplus = w5.plus(u);				   // <--------------
  cout << w5 << ".plus" << u << " = " << wXplus << endl;

// ***** Various norms and comparisons

	// The invariant square of difference between two 4-vectors
	//--------------------------------------------------------

  d = w1.diff2(w2);					   // <--------------
  cout << w1 << ".diff2" << w2 << " = " << d << endl;

	// Although this is Lorentz invariant, its sign does
	// depend on the metric.  The euclidean differnece squared
	// does not:

  d = w1.delta2Euclidean(w2);				   // <--------------
  cout << w1 << ".delta2Euclidean" << w2
				<< " = " << d << endl;

	// Here is how to test relative nearness in Euclidean sense
	//---------------------------------------------------------

  w5 = w1 + .000001 * w3;

  d = w5.howNear(w1);					   // <--------------
  cout << w5 << ".howNear" << w1 << " = " << d << endl;

	// This is the euclidean norm of the difference, made scale
	// invariant by dividing by sqrt (v1.dot(v2) + average_t**2)

  d = (100*w5).howNear(100*w1);
  cout << 100*w5 << ".howNear" << 100*w1 << " = " << d << endl;

  if ( ! w5.isNear(w1) ) {				   // <--------------
    cout << w5 << " is not near " << w1 << endl;
  }
	// isNear compares to the tolerance, which has a
	// default value of 100 times Float8 machine epsilon:
	// about 2.2E-14.

  if ( w5.isNear(w1, .01) ) {				   // <--------------
    cout << "With tolerance .01, " << w5
		<< " is near " << w1 << endl;
  }

	// You can control the tolerance in each instance or
	// change the default tolerance

  Float8 oldTol = LorentzVector::setTolerance (.001);	   // <--------------

  if ( w5.isNear(w1) ) {
    cout << "With tolerance " << w3.getTolerance()
	<< ",\n" <<  w5 << " is near " << w1 << endl;
  }

  LorentzVector::setTolerance (oldTol);
	// I leave things the way I found them

	// Here is how to test relative nearness in Center of Mass frame
	//--------------------------------------------------------------

	// The above definition is not Lorentz invariant, but may
	// sensibly be applied to any pair of LorentzVectors.  An
	// intuitively attractive alternative, for timelike vectors,
	// is the Lorentz invariant version:

  w5 = w1 + .01 * w3;

  d = w5.howNearCM(w1);					   // <--------------

  	// This boosts both vectors into their joint CM frame
	// and applies howNEar.  Thus howNearCM() is both scale
	// invariant and Lorentz invariant.

  cout << w5 << ".howNearCM" << w1 << " = " << w5.howNearCM(w1) << endl;

  w6 = w5 * 1000;
  w7 = w1 * 1000;
  cout << w6 << ".howNearCM" << w7 << " = " << w6.howNearCM(w7) << endl;

  w6 = boostOf(w5, Z_HAT, .995);
  w7 = boostOf(w1, Z_HAT, .995);
  cout << w6 << ".howNearCM" << w7 << " = " << w6.howNearCM(w7) << endl;

	// The same options for tolerance are available for
	// isNearCM() as are available for isNear()

  if ( w5.isNearCM(w1, .03) ) {				   // <--------------
    cout << "With tolerance .03, " << w5
		<< " is near " << w1 << endl;
  }

  return 0;

} /* end of main */