// lvtest.cc -- general test of all features of LorentzVector class
//-----------------------------------------------------------------

// Modifications since 2/1/98:
// ---------------------------
//
//  2/ 6/98 Dictionary ordering is z first:  Changed test vectors in comparisons
//	    in dictionary ordering section.
//  2/19/98 Modifications to break up into a dozen smaller functions for
//	    compilation speed.
//  3/23/98 Inclusion of tests of quantitative nearness methods
//
//  5/27/98 Replace tests for exact equality with eq which tolerates 10**-15
//	    to accomodate x86 math inaccuracy
//  6/15/98 Added namespace support
//

#include "ZMutility/ZMenvironment.h"

#include "PhysicsVectors/SpaceVector.h"
#include "PhysicsVectors/UnitVector.h"
#include "PhysicsVectors/LorentzVector.h"
ZM_USING_NAMESPACE( zmpv )

#include "ZMutility/cmath"
using CMATH_NAMESPACE::fabs;

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

#include <vector>
using std::vector;

#include "ZMutility/PI.h"
const Float8 PI =  ZMpi;


// Variables used by all the parts:

bool problems (false);
LorentzVector vans;
LorentzVector v1 ( 1., 2., 3., Tcomponent(4.) );
LorentzVector v2 ( Tcomponent(5), 1., 2, 3. );
const LorentzVector v1c ( Tcomponent(6), 2, 3., 4. );
SpaceVector vspace (5,7,2);
LorentzVector v3( vspace, 30);
LorentzVector v4( 20, vspace);
LorentzVector v5( vspace, Tcomponent(40) );
LorentzVector v6( Tcomponent(50), vspace );
LorentzVector v7( Tcomponent(25) );
LorentzVector v8( -10 );
SpaceVector vg = v1.v();
SpaceVector svans(1,2,3);
SpaceVector vx ( 3, -4, 12);
const LorentzVector v3c( 3, 5, 7, Tcomponent(9) );
const LorentzVector v4c( 6, 2, 8, Tcomponent(7) );
double toler;

// The parts are:

void constructors();
void accessors();
void compAccessors();
void modifiers();
void setMultiple();
void arithmetic();
void comparisons();
void tolerance();
void nearCM();
void nearness();
void combining();
void intrinsic();
void kinematic();
void rotations();
void rotations2();
void rotations3();
void boosts();
void units();
void input();



bool checkAxialRotation ( const LorentzVector &x,
			  const SpaceVector &d,
			  const Float8 &phi,
			  const LorentzVector& result );


bool close ( double a, double b ) {
  return ( fabs(a-b) < 1.0E-10 );
}
bool eq ( double a, double b ) {
  return ( fabs(a-b) < 1.0E-15 );
}

double fNear( LorentzVector w1, LorentzVector w2 ) {
  return sqrt( w1.delta2Euclidean(w2)/(w1.euclideanNorm()*w2.euclideanNorm()));
}



int main() {

// For purposes of this test, it is safe and necessary to ignore ZMexceptions.
// This test exercises some of them but we want to continue beyond each one.
  ZMxPhysicsVectors::setHandler( ZMexIgnoreAlways() );


  cout << "\n\n";
  cout << "   *************\n";
  cout << "   LorentzVector\n";
  cout << "   *************\n";

  constructors();
  accessors();
  compAccessors();
  modifiers();
  setMultiple();
  comparisons();
  tolerance();
  nearCM();
  nearness();
  combining();
  intrinsic();
  kinematic();
  rotations();
  rotations2();
  rotations3();
  boosts();
  units();
  input();

  if (problems) {
    cout << "\n ?????????????????? -- End of tests:  PROBLEMS DETECTED\n\n";
    return -1;
  } else {
    cout << "\nEnd of tests -- All automated checking passed.\n\n";
    return 0;
  }

} /* End of lvtest main() */


void constructors()  {

  cout << "\n --- Constructors and Output\n\n";

	cout << "{ Cartesian Constructor }\n";

#ifdef ZAP_timeComp
  LorentzVector v1 ( 1., 2., 3., 4. );
// Properly Fails to compile:
// call of overloaded constructor
// `LorentzVector(double, double, double, double)' is ambiguous
//  candidates are: LorentzVector(Tcomponent, double, double, double)
// 		    LorentzVector(double, double, double, Tcomponent)
#endif // ZAP_timeComp

  if ( (v1.x() == 1.) && (v1.y()==2) && 3==v1.z() && v1.t()==4. ) {
    cout << "OK: constructed " << v1 << endl;
  } else {
    cout << "????? constructed " << v1 << endl;
    problems=true;
  }

  if ( (v2.x() == 1.) && (v2.y()==2) && 3==v2.z() && v2.t()==5 ) {
    cout << "OK: constructed using int Tcomponent " << v2 << endl;
  } else {
    cout << "????? constructed " << v2 << endl;
    problems=true;
  }

  if ( sizeof(v1) == 32 ) {
    cout << "OK: sizeof LorentzVector is " << sizeof (v1) << endl;
  } else {
    cout << "????? sizeof LorentzVector is " << sizeof (v1) << endl;
    problems=true;
  }

  if ( (v1c.z() == 4.) && (v1c.x()==2) && 3==v1c.y() && v1c.t()==6 ) {
    cout << "OK: constructed const " << v1c << endl;
  } else {
    cout << "????? constructed const " << v1c << endl;
    problems=true;
  }
	cout << "{ V,t Constructor }\n";

  if ( (v3.x() == 5.) && (v3.y()==7) && 2==v3.z() && v3.t()==30 ) {
    cout << "OK: constructed using v,t " << v3 << endl;
  } else {
    cout << "????? constructed " << v3 << endl;
    problems=true;
  }

  if ( (v4.x() == 5.) && (v4.y()==7) && 2==v4.z() && v4.t()==20 ) {
    cout << "OK: constructed using t, v " << v4 << endl;
  } else {
    cout << "????? constructed " << v4 << endl;
    problems=true;
  }

  if ( (v5.x() == 5.) && (v5.y()==7) && 2==v5.z() && v5.t()==40 ) {
    cout << "OK: constructed using v, Tcomponent(t) " << v5 << endl;
  } else {
    cout << "????? constructed " << v5 << endl;
    problems=true;
  }

  if ( (v6.x() == 5.) && (v6.y()==7) && 2==v6.z() && v6.t()==50 ) {
    cout << "OK: constructed using Tcomponent(t), v " << v6 << endl;
  } else {
    cout << "????? constructed " << v6 << endl;
    problems=true;
  }

	cout << "{ t-only Constructor }\n";

  SpaceVector temp(v7.v());
  if ( temp == SpaceVector(0,0,0) && v7.t()==25 ) {
    cout << "OK: constructed using Tcomponent(int) " << v7 << endl;
  } else {
    cout << "????? constructed " << v7 << endl;
    problems=true;
  }

  if ( v8.getV() == SpaceVector(0,0,0) && v8.t()== -10. ) {
    cout << "OK: constructed using t " << v8 << endl;
  } else {
    cout << "????? constructed " << v8 << endl;
    problems=true;
  }

} // end of constructors();


void accessors()  {

// Accessors

  cout << "\n --- Accessors\n\n";

  // Straight x y z

	cout << "{ get Cartesian components }\n";
  double xx = v1.x();
  double yy = v1.y();
  double zz = v1.z();
  double tt = v1.t();
  if ( (xx==1) && (yy==2) && (zz==3) && (tt==4) ) {
    cout << "OK: v1 = " << v1 << endl;
  } else {
    cout << "????? v1 = " << v1 << endl;
    problems=true;
  }

  xx = v1c.x();
  yy = v1c.y();
  zz = v1c.z();
  tt = v1c.t();
  if ( (xx==2) && (yy==3) && (zz==4) && (tt==6) ) {
    cout << "OK: v1c = " << v1c << endl;
  } else {
    cout << "????? v1c = " << v1c << endl;
    problems=true;
  }

	cout << "{ output Cartesian components }\n";

  cout << "    ... ostream<< x, y, z, t ==> "
	<< v1.x() << ", " << v1.y() << ", " << v1.z() << ", " << v1.t() << endl;

  cout << "    ... ostream<< x, y, z, t  ==> "
	<< v1c.x() << ", " << v1c.y() << ", " << v1c.z() << ", " << v1c.t() << endl;

	cout << "{ use Cartesian components in expressions }\n";

  // NOTE: When the constructor of LorentzVector(Scalar) was not explicit,
  // gcc misbehaved in this case, converting 6 to a LV, multiplying by the
  // double formed from v1.z(), comverting 18 to a LV, and applying the
  // LV::operator==(LV).  Explicit cures this.
  if ( 6*v1.z() == 18 ) {
    cout << "OK: 6*v1.z() = " << 6*v1.z() << endl;
  } else {
    cout << "????? 6*v1.z() = " << 6*v1.z() << endl;
    problems=true;
  }

  if ( 6+v1.z() == 9 ) {
    cout << "OK: 6+v1.z() = " << 6+v1.z() << endl;
  } else {
    cout << "????? 6+v1.z() = " << 6+v1.z() << endl;
    problems=true;
  }

  if ( v1.z()*7 == 21 ) {
    cout << "OK: v1.z()*7 = " << v1.z()*7 << endl;
  } else {
    cout << "????? v1.z()*7 = " << v1.z()*7 << endl;
    problems=true;
  }

  if ( 6*v1.z() == 18 ) {
    cout << "OK: 6*v1.z() = " << v1.z()*18 << endl;
  } else {
    cout << "????? 6*v1.z() = " << 6*v1.z() << endl;
    problems=true;
  }

  if ( v1.t()+3 == 7 ) {
    cout << "OK: v1.t()+3 = " << v1.t()+3 << endl;
  } else {
    cout << "????? v1.t()+3 = " << v1.t()+3 << endl;
    problems=true;
  }

  if ( 5*v1c.y() == 15 ) {
    cout << "OK: 5*v1c.z() = " << 5*v1c.z() << endl;
  } else {
    cout << "????? 5*v1c.z() = " << 5*v1c.z() << endl;
    problems=true;
  }

  if ( v1c.y()/2 == 1.5 ) {
    cout << "OK: v1c.y()/2 = " << v1c.y()/2 << endl;
  } else {
    cout << "????? v1c.y()/2 = " << v1c.y()/2 << endl;
    problems=true;
  }

  // Observation:  v.x_ directly does not compile, as desired.
#ifdef ZAP1
  cout << v.x_;
#endif // ZAP1

} // end of accessors

void compAccessors() {

  // Comparison of accessors

	cout << "{ compare Cartesian components }\n";
  if ( v1.x() == 1  ) {
    cout << "OK: Equality comparison on an accessor works\n";
  } else {
    cout << "????? Equality comparison on an accessor FAILS!!!!!\n";
    problems=true;
  }

  if ( v1.x() != v1.y() ) {
    cout << "OK: comp != comp works\n";
  } else {
    cout << "????? comp != comp FAILS!!!!\n";
    problems=true;
  }

  if ( v1.x() == v1.y() ) {
    cout << "????? comp == comp FAILS!!!!\n";
    problems=true;
  } else {
    cout << "OK: comp == comp works\n";
  }

  if ( v1.x() > v1.y() ) {
    cout << "?????comp > comp FAILS!!!!\n";
    problems=true;
  } else {
    cout << "OK: comp > comp works\n";
  }


  if ( v1.y() < v1.x() ) {
    cout << "????? comp < comp FAILS!!!!\n";
    problems=true;
  } else {
    cout << "OK: comp < comp works\n";
  }


  double xeq = 4.0;

  if ( xeq == v1.z() ) {
    cout << "????? Scalar == comp FAILS!!!!\n";
    problems=true;
  } else {
    cout << "OK: Scalar == comp works\n";
  }

  const double ceq = 4.0;
  if ( ceq == v1.z() ) {
    cout << "????? Scalar == comp FAILS!!!!\n";
    problems=true;
  } else {
    cout << "OK: Scalar == comp works\n";
  }

  if ( v1c.x() == 2 ) {
    cout << "OK: Equality comparison on a const vector accessor works\n";
  } else {
    cout << "?????Equality comparison on a const vector accessor FAILS!!!!!\n";
    problems=true;
  }

  if ( v1.x() != v1.y() ) {
    cout << "OK: comp != comp works\n";
  } else {
    cout << "?????comp != comp FAILS!!!!\n";
    problems=true;
  }

  if ( v1.x() == v1.y() ) {
    cout << "????? comp == comp FAILS!!!!\n";
    problems=true;
  } else {
    cout << "OK: comp == comp works\n";
  }

  if ( v1.x() > v1.y() ) {
    cout << "????? comp > comp FAILS!!!!\n";
    problems=true;
  } else {
    cout << "OK: comp > comp works\n";
  }

  if ( v1.y() < v1.x() ) {
    cout << "????? comp < comp FAILS!!!!\n";
    problems=true;
  } else {
    cout << "OK: comp < comp works\n";
  }

  xeq = 4.0;


  if ( v1c.y() < v1c.z() ) {
    cout << "OK: const comp == const comp works\n";
  } else {
    cout << "????? const comp = const comp FAILS!!!!\n";
    problems=true;
  }



	cout << "{ Compare v component }\n";

  if ( vg == svans ) {
    cout << "OK: v1.v() = " << v1.v() << endl;
  } else {
    cout << "????? v1.v() = " << v1.v() << endl;
    problems=true;
  }

  if ( svans == v1.v() ) {
    cout << "OK: direct lhs use of v1.v() in comparison works " << endl;
  } else {
    cout << "????? direct lhs use of v1.v() in comparison FAILS " << endl;
    problems=true;
  }

  if ( v1.v() == svans ) {
    cout << "OK: direct rhs use of v1.v() in comparison works " << endl;
  } else {
    cout << "????? direct rhs use of v1.v() in comparison FAILS " << endl;
    problems=true;
  }

  svans.set(1,2,2);

  if ( svans != v1.v() ) {
    cout << "OK: direct lhs use of v1.v() in != comparison works " << endl;
  } else {
    cout << "????? direct lhs use of v1.v() in != comparison FAILS " << endl;
    problems=true;
  }

  if ( v1.v() != svans ) {
    cout << "OK: direct rhs use of v1.v() in != comparison works " << endl;
  } else {
    cout << "????? direct rhs use of v1.v() in != comparison FAILS " << endl;
    problems=true;
  }

  v7 = v1;

  if ( v1.v() == v7.v() ) {
    cout << "OK: direct both-side use of .v() in == comparison works "
				<< endl;
  } else {
    cout << "????? direct both-side use of v1.v() in == comparison FAILS "
				<< endl;
    problems=true;
  }

  // Note that direct comparison of accessors for V using dictionary ordering
  // is not supported:  Nobody will have a list of ACCESSORS to sort.

  // Setting different components

	cout << "{ Set Cartesian components }\n";

  v6.setX( 24 );
  v6.setY( 3. );
  v6.setZ( 100 );
  v6.setT( 55 );
  vans.set(24,3,100,Tcomponent(55));
  if (v6==vans) {
    cout << "OK: set x=24, y=3, z=100, t=55 ==> " << v6 << endl;
  } else {
    cout << "????? set x=24, y=3, z=100, t=55 ==> " << v6 << endl;
    problems=true;
  }

	cout << "{ Set SpaceVector component }\n";

  v5.setT( 90 );
  v5.setV( 2*v6.v() );
  vans.set(48,6,200,Tcomponent(90));
  if (v5==vans) {
    cout << "OK: set v=2*(24,3,100), t=90 ==> " << v5 << endl;
  } else {
    cout << "????? set v=2*(24,3,100), t=90 ==> " << v5 << endl;
    problems=true;
  }


#ifdef ZAP2
  v1c.v() = .05;
  // Properly fails, with message
  //		`WontAlter_V_of_a_const' does not define operator=
#endif

} // end of compAccessors()


void modifiers()  {

  cout << "\n --- Component Modify/Assigners\n\n";

// Component modifiers ( v.x() += Scalar, etc.)

	cout << "{ Modifying cartesian components }\n";

  v1.set(1,2,3,Tcomponent(4));

  cout << v1 << " Add 5 to X component\n";
  v1.setX( v1.x() + 5 );
  vans.set(6,2,3,Tcomponent(4));
  if (v1==(vans)) {
    cout << "OK: " << v1 << " ==> " << vans << endl;
  } else {
    cout << "????? " << v1 << " ==> " << vans << endl;
    problems = true;
  }

  cout << v1 << " v1.y() *= v1.x\n";
  v1.setY( v1.y() * v1.x() );
  vans.set(6,12,3,Tcomponent(4));
  if (v1==(vans)) {
    cout << "OK: " << v1 << " ==> " << vans << endl;
  } else {
    cout << "????? " << v1 << " ==> " << vans << endl;
    problems = true;
  }

  cout << v1 << " Subtract 5 from Z component\n";
  v1.setZ( v1.z() - 5 );
  vans.set(6,12,-2,Tcomponent(4));
  if (v1==(vans)) {
    cout << "OK: " << v1 << " ==> " << vans << endl;
  } else {
    cout << "????? " << v1 << " ==> " << vans << endl;
    problems = true;
  }

  cout << v1 << " v1.t() /= 2\n";
  v1.setT( v1.t() / 2 );
  vans.set(6,12,-2,Tcomponent(2));
  if (v1==(vans)) {
    cout << "OK: " << v1 << " ==> " << vans << endl;
  } else {
    cout << "????? " << v1 << " ==> " << vans << endl;
    problems = true;
  }

	cout << "{ Modifying V component }\n";

  cout << v1 << " v1.v() /= 2\n";
  v1.setV( v1.v()/2 );
  vans.set(3,6,-1,Tcomponent(2));
  if (v1==(vans)) {
    cout << "OK: " << v1 << " ==> " << vans << endl;
  } else {
    cout << "????? " << v1 << " ==> " << vans << endl;
    problems = true;
  }

  cout << v1 << " v1.v() += (3,4,5)\n";
  v1.setV( v1.v()+SpaceVector(3,4,5) );
  vans.set(6,10,4,Tcomponent(2));
  if (v1==(vans)) {
    cout << "OK: " << v1 << " ==> " << vans << endl;
  } else {
    cout << "????? " << v1 << " ==> " << vans << endl;
    problems = true;
  }

} // end of modifiers


void setMultiple() {

  cout << "\n --- Setting multiple components at once\n\n";

	cout << "{ Cartesian components }\n";

  v1.set ( 5, 7, 9, Tcomponent(11));
  if ( (v1.x() == 5) && (v1.y() == 7) && (v1.z() == 9) && v1.t() == 11) {
    cout << "OK: " << v1 << endl;
  } else {
    cout << "????? " << v1 << endl;
    problems = true;
  }

  v1.set ( Tcomponent(12), 15, 17, 19);
  if ( (v1.x() == 15) && (v1.y() == 17) && (v1.z() == 19) && v1.t() == 12) {
    cout << "OK: " << v1 << endl;
  } else {
    cout << "????? " << v1 << endl;
    problems = true;
  }

	cout << "{ V, t components }\n";

  svans.set (3, 4, 5);
  v1.set ( Tcomponent(10), svans );
  if ( (v1.x() == 3) && (v1.y() == 4) && (v1.z() == 5) && v1.t() == 10) {
    cout << "OK: " << v1 << endl;
  } else {
    cout << "????? " << v1 << endl;
    problems = true;
  }

  svans.set (3, -4, 5);
  v1.set ( svans, 8 );
  if ( (v1.x() == 3) && (v1.y() == -4) && (v1.z() == 5) && v1.t() == 8) {
    cout << "OK: " << v1 << endl;
  } else {
    cout << "????? " << v1 << endl;
    problems = true;
  }

} // end of setMultiple()

void arithmetic() {

// Arithmetic operations

  cout << "\n --- Arithmetic\n\n";

	cout << "{ + and - }\n";

  v1.set( 3, 5, 7, Tcomponent(8) );
  v2.set( 6, 2, 8, Tcomponent(4));

  v3 = v1 + v2;
  vans.set ( 9, 7, 15, Tcomponent(12) );
  if (v3.isNear(vans)) {
    cout << "OK: " << v1 << " + " << v2 << " = " << v3 << endl;
  } else {
    cout << "????? " << v1 << " + " << v2 << " = " << v3 << endl;
    problems = true;
  }

  v3 = v1 - v2;
  vans.set ( -3, 3, -1, Tcomponent(4) );
  if (v3.isNear(vans)) {
    cout << "OK: " << v1 << " - " << v2 << " = " << v3 << endl;
  } else {
    cout << "????? " << v1 << " - " << v2 << " = " << v3 << endl;
    problems = true;
  }

  v3 = - v2;
  vans.set ( -6, -2, -8, Tcomponent(-4) );
  if (v3.isNear(vans)) {
    cout << "OK: " << "- " << v2 << " = " << v3 << endl;
  } else {
    cout << "????? " << "- " << v2 << " = " << v3 << endl;
    problems = true;
  }

  v3 = v4c + v3c;
  vans.set ( 9, 7, 15, Tcomponent(16) );
  if (v3.isNear(vans)) {
    cout << "OK: " << v4c << " + " << v3c << " = " << v3 << endl;
  } else {
    cout << "????? " << v4c << " + " << v3c << " = " << v3 << endl;
    problems = true;
  }

  v3 = v4c - v3c;
  vans.set ( 3, -3, 1, Tcomponent(-2) );
  if (v3.isNear(vans)) {
    cout << "OK: " << v4c << " - " << v3c << " = " << v3 << endl;
  } else {
    cout << "????? " << v4c << " - " << v3c << " = " << v3 << endl;
    problems = true;
  }

  v3 = - v3c;
  vans.set ( -3, -5, -7, Tcomponent(-9) );
  if (v3.isNear(vans)) {
    cout << "OK: " << "- " << v3c << " = " << v3 << endl;
  } else {
    cout << "????? " << "- " << v3c << " = " << v3 << endl;
    problems = true;
  }

	cout << "{ * and / }\n";

  double pointfive = .5;
  v3 = v2 * pointfive;
  vans.set ( 3, 1, 4, Tcomponent(2) );
  if (v3==vans) {
    cout << "OK: " << v2 << " * .5 = " << v3 << endl;
  } else {
    cout << "????? " << v2 << " * .5 = " << v3 << endl;
    problems = true;
  }

  v3 = v2 / 4;
  vans.set ( 1.5, .5, 2, Tcomponent(1.0) );
  if (v3==vans) {
    cout << "OK: " << v2 << " / 4 = " << v3 << endl;
  } else {
    cout << "????? " << v2 << " / 4 = " << v3 << endl;
    problems = true;
  }

  v3 = 1.5 * v2;
  vans.set ( Tcomponent(6), 9, 3, 12 );
  if (v3==vans) {
    cout << "OK: 1.5 * " << v2 << " = " << v3 << endl;
  } else {
    cout << "????? 1.5 * " << v2 << " = " << v3 << endl;
    problems = true;
  }

  v3 = v4c * .5;
  vans.set ( Tcomponent(3.5), 3, 1, 4 );
  if (v3==vans) {
    cout << "OK: " << v4c << " * .5 = " << v3 << endl;
  } else {
    cout << "????? " << v4c << " * .5 = " << v3 << endl;
    problems = true;
  }

  v3 = v4c / .5;
  vans.set ( 12, 4, 16, Tcomponent(14) );
  if (v3==vans) {
    cout << "OK: " << v4c << " / .5 = " << v3 << endl;
  } else {
    cout << "????? " << v4c << " / .5 = " << v3 << endl;
    problems = true;
  }

  v3 = 2.5 * v4c;
  vans.set (15, 5, 20, Tcomponent(17.5) );
  if (v3==vans) {
    cout << "OK: 2.5 * " << v4c << " = " << v3 << endl;
  } else {
    cout << "????? 2.5 * " << v4c << " = " << v3 << endl;
    problems = true;
  }

	cout << "{ Linear expression }\n";

  v3.set(9,3,12, Tcomponent(13.5));
  v4 = v3 + 5 * v1 + v4c/2;
  vans.set (27, 29, 51, Tcomponent(57));
  if (v4==vans) {
    cout << "OK: " << v3 << " + 5 * " << v1 << endl <<
	"            + " << v4c << " / 2 = " <<  v4 << endl;
  } else {
    cout << "????? " << v3 << " + 5 * " << v1 << endl <<
	"            + " << v4c << " / 2 = " <<  v4 << endl;
    problems = true;
  }

	cout << "{ Modify/assign on whole vectors }\n";

  v1.set ( 4, 8, 6, Tcomponent(-10) );
  cout << "    ... v1 = " << v1 << endl;
  v1 *= .5;
  vans.set ( 2, 4, 3, Tcomponent(-5) );
  if (v1==vans) {
    cout << "OK: v1 *= .5 ==> " << v1 << endl;
  } else {
    cout << "????? v1 *= .5 ==> " << v1 << endl;
    problems = true;
  }

  v1 /= .1;
  vans.set ( 20, 40, 30, Tcomponent(-50) );
  if (v1.isNear(vans)) {
    cout << "OK: v1 /= .1 ==> " << v1 << endl;
  } else {
    cout << "????? v1 /= .1 ==> " << v1 << endl;
    problems = true;
  }

  v1 += LorentzVector ( 5, -5, 10, Tcomponent(20) );
  vans.set ( 25, 35, 40, Tcomponent(-30) );
  if (v1.isNear(vans)) {
    cout << "OK: v1 += (5, -5, 10) ==> " << v1 << endl;
  } else {
    cout << "????? v1 += (5, -5, 10) ==> " << v1 << endl;
    problems = true;
  }

  v1 -= v1/5;
  vans.set ( 20, 28, 32, Tcomponent(-24) );
  if (v1.isNear(vans)) {
    cout << "OK: v1 -= v1/5 ==> " << v1 << endl;
  } else {
    cout << "????? v1 -= v1/5 ==> " << v1 << endl;
    problems = true;
  }

} // end of arithmetic()

void comparisons() {

// Comparisons

  cout << "\n --- Comparisons\n\n";

  v1.set (2, 4,  6, Tcomponent(8));
  v2.set (1, 5, 9, Tcomponent(-2));
  // Previously:
  //  const LorentzVector v3c( 3, 5, 7, 9);
  //  const LorentzVector v4c( 6, 2, 8, 7);
  v5 = v1;

	cout << "{ == and != }\n";

  if ( ! (v1 == v5) ) {
    cout << "????? v1 == v5 false  erroneously! v1 = "
	<< v1 << " v5 = " << v5 << endl;
    problems = true;
  }
  cout << "OK: " << v1 << " == " << v5 << endl;

  if ( v1 == v2 ) {
    cout << "????? v1 == v2 true erroneously! v1 = "
	<< v1 << " v2 = " << v2 << endl;
    problems = true;
  }
  cout << "OK: " << "NOT " << v1 << " == " << v2 << endl;

  if ( v1 != v5 ) {
    cout << "????? v1 != v5 true erroneously! v1 = "
	<< v1 << " v5 = " << v5 << endl;
    problems = true;
  }
  cout << "OK: " << "NOT " << v1 << " != " << v5 << endl;

  if ( ! (v1 != v2) ) {
    cout << "????? v1 != v2 false erroneously! v1 = "
	<< v1 << " v2 = " << v2 << endl;
    problems = true;
  }
  cout << "OK: " << v1 << " != " << v2 << endl;

  if ( v1 == v3c ) {
    cout << "????? v1 == v3c true erroneously! v1 = "
	<< v1 << " v3c = " << v3c << endl;
    problems = true;
  }
  cout << "OK: " << "NOT " << v1 << " == " << v2 << endl;

  if ( ! (v3c != v4c) ) {
    cout << "????? v3c != v4c false erroneously! v3c = "
	<< v3c << " v4c = " << v4c << endl;
    problems = true;
  }
  cout << "OK: " << v3c << " != " << v4c << endl;


	cout << "{ Dictionary Ordering: >, <, >=, <= }\n";

  v1.set (6, 4, 2, Tcomponent(8) );	// CHANGED 2/6/98 FROM (2,4,6,8)
  v2.set (9, 5, 3, Tcomponent(7));	// CHANGED 2/6/98 FROM (3,5,9,7)
  v5.set (3, 5, 2, Tcomponent(8) );	// CHANGED 2/6/98 FROM (2,5,3,8)
  v6 = v1;

  // Previously:
  //  const LorentzVector v3c( 3, 5, 7, 9);
  const LorentzVector v33c( 7, 5, 3, Tcomponent(9)); // CHANGED 2/6/98


  if ( ! (v1 <  v5) ) {
    cout << "????? v1 < v5 false erroneously! v1 = "
	<< v1 << " v5 = " << v5 << endl;
    problems = true;
  }
  cout << "OK: " << v1 << " <  " << v5 << endl;

  if ( v1 <  v6 ) {
    cout << "????? v1 < v6 true erroneously! v1 = "
	<< v1 << " v6 = " << v6 << endl;
    problems = true;
  }
  cout << "OK: " << "NOT " << v1 << " <  " << v6 << endl;

  if ( v1 >  v6 ) {
    cout << "????? v1 > v6 true erroneously! v1 = "
	<< v1 << " v6 = " << v6 << endl;
    problems = true;
  }
  cout << "OK: " << "NOT " << v1 << " >  " << v6 << endl;

  if ( v5 >=  v3c ) {
    cout << "????? v5 >= v3c true erroneously! v5 = "
	<< v5 << " v5 = " << v3c << endl;
    problems = true;
  }
  cout << "OK: " << "NOT " << v5 << " >=  " << v3c << endl;

   if ( ! (v1 <= v5 ) ) {
    cout << "????? v1 <= v5 false erroneously! v1 = "
	<< v1 << " v5 = " << v5 << endl;
    problems = true;
  }
  cout << "OK: " << v1 << " <= " << v5 << endl;

   if ( ! (v1 <= v6 ) ) {
    cout << "????? v1 <= v6 false erroneously! v1 = "
	<< v1 << " v6 = " << v6 << endl;
    problems = true;
  }
  cout << "OK: " << v1 << " <= " << v6 << endl;

   if ( ! (v6 <= v1 ) ) {
    cout << "????? v6 <= v1 true erroneously! v6 = "
	<< v6 << " v1 = " << v1 << endl;
    problems = true;
  }
  cout << "OK: " << v6 << " <= " << v1 << endl;

} // end of comparisons();


void tolerance() {

// Comparisons within tolerance

  cout << "\n --- IsNear() and isNearCM(): Comparisons within tolerance\n\n";

	cout << " { isNear() with default tolerance }\n";

  toler = LorentzVector::getTolerance();
  if ( (toler > 2.e-14) && (toler < 3.e-14) ) {
    cout << "OK: Default tolerance starts at (expected 2.22e-14) ==> "
							<< toler << endl;
  } else {
    cout << "????? Default tolerance starts at (expected 2.22e-14) ==> "
							<< toler << endl;
    problems=true;
  }
  LorentzVector::setTolerance(toler);

  v3 = v1 + .2 * toler * v2;
  if (v3 != v1) {
    cout << "OK: " << "adding .2 * tolerance * times " << v2 << endl <<
	"   to v1 = " << v1 << " gives V3 != V1\n";
  } else {
    cout << "????? " << "adding .2 * tolerance * times " << v2 << endl <<
	"   to v1 = " << v1 << " gives V3 != V1\n";
    problems = true;
  }

  if ( v3 .isNear(v1) ) {
    cout << "OK: " << "Added .20 * tolerance * v2 -- isNear\n";
  } else {
    cout << "????? Added .2 * tolerance -- isNear FALSE\n";
    problems = true;
  }

  v3 = v1 + .7 * toler * v2;
  if ( v1 .isNear(v3) ) {
    cout << "OK: Added .72 * tolerance * v2 -- isNear\n";
  } else {
    cout << "????? " << "Added .72 * tolerance * v2 -- isNear FALSE\n";
    problems = true;
  }

  v3 = v1 + .85 * toler * v2;
  if ( v1 .isNear(v3) ) {
    cout << "OK: Added .85 * tolerance * v2 -- isNear\n";
  } else {
    cout << "????? " << "Added .85 * tolerance * v2 -- isNear FALSE\n";
    problems = true;
  }

  v3 = v1 + .86 * toler * v2;
  if ( v1 .isNear(v3) ) {
    cout << "????? " << "Added .86 * tolerance * v2 -- isNear\n";
    problems = true;
  } else {
    cout << "OK: Added .86 * tolerance * v2 -- isNear FALSE\n";
  }

  v3 = v1 + 1.2 * toler * v2;
  if (  v3.isNear(v1) ) {
    cout << "????? " << "Added 1.2 * tolerance * v2 -- isNear\n";
    problems = true;
  } else {
    cout << "OK: Added 1.2 * tolerance * v2 -- isNear FALSE\n";
  }

  v3 = v1 + 1.8 * toler * v2;
  if (  v3.isNear(v1) ) {
    cout << "????? Added 1.8 * tolerance * v2 -- isNear TRUE\n";
    problems = true;
  } else {
    cout << "OK: " << "Added 1.2 * tolerance * v2 -- isNear FALSE\n";
  }


	cout << " { isNear() with specified tolerance }\n";


  v1.set ( 1, 3, 6, Tcomponent(5) );
  double m = sqrt(v1.getV().mag2() + v1.t()*v1.t() );
  cout << "    .... dealing with a vector " << v1
				<< " of magnitude " << m << endl;
  double small = .01;
  v2 = v1 + small * LorentzVector( 1, -1, 1, Tcomponent(-1) );

  if ( v1.isNear(v2) ) {
    cout << "?????  with (default) tolerance " << toler << ", "
	<< v1 << " isNear " << v2 << endl;
    problems = true;
  } else {
    cout << "OK: " << "with (default) tolerance " << toler << ",\n"
	"        " << v1 << " is not Near " << v2 << endl;
  }

  if ( v1.isNear(v2 ,  .5*small/m) ) {
    cout << "????? with tolerance " << .5*small/m <<
	", " << v1 << " isNear " << v2 << endl;
    problems = true;
  } else {
    cout << "OK: " << "with tolerance " << .5*small/m << ", "
	<< v1 << "\n                is not Near " << v2 << endl;
  }

  if ( v1.isNear(v2 , 1.5*small/m) ) {
    cout << "?????with tolerance " << 1.5*small/m << ", "
	<< v1 << " isNear " << v2 << endl;
    problems = true;
  } else {
    cout << "OK:  with tolerance " << 1.5*small/m << ", "
    	<< v1 << "\n                is not Near " << v2 << endl;
  }

  if ( v1.isNear(v2 , 2.5*small/m) ) {
    cout << "OK:  with tolerance " << 2.5*small/m << ", "
	<< v1 << "\n                isNear " << v2 << endl;
  } else {
    cout << "????? with tolerance " << 2.5*small/m << ", "
    	<< v1 << " is not Near " << v2 << endl;
    problems = true;
  }


  if ( v1.isNear(v2 , 3.5*small/m) ) {
    cout << "OK:  " << "with tolerance " << 3.5*small/m << ", "
	<< v1 << "\n                isNear " << v2 << endl;
  } else {
    cout << "?????with tolerance " << 3.5*small/m << ", "
    	<< v1 << " is not Near " << v2 << endl;
    problems = true;
  }

	cout << " { setTolerance }\n";

  LorentzVector::setTolerance ( .1 );
  cout << "    ... tolerance set to .1\n";

  LorentzVector wj (1, 1, -1, Tcomponent(-1));

  v2 = v1 + .02 * wj;

  if ( v1.isNear(v2) ) {
    cout << "OK: " <<
		"with new default tolerance " << v1 <<
		"\n                isNear " << v2 << endl;
  } else {
    cout << "?????with new default tolerance "
	<< v1 << " is not Near " << v2 << endl;
    problems = true;
  }

  v2 = v1 + .05 * wj;

  if ( v1.isNear(v2) ) {
    cout << "OK: " <<
	"with new default tolerance " << v1 <<
	"\n                isNear " << v2 << endl;
  } else {
    cout << "?????with new default tolerance "
	<< v1 << " is not Near " << v2 << endl;
    problems = true;
  }

  v2 = v1 + .1 * wj;

  if ( v1.isNear(v2) ) {
    cout << "OK: with new default tolerance " << v1 <<
	"\n                isNear " << v2 << endl;
  } else {
    cout << "????? with new default tolerance " <<
			v1 << " is not Near " << v2 << endl;
    problems = true;
  }

  v2 = v1 + .2 * wj;

  if ( v1.isNear(v2) ) {
    cout << "OK: with new default tolerance " << v1 <<
	"\n                isNear " << v2 << endl;
  } else {
    cout << "?????with new default tolerance " <<
	v1 << " is not Near " << v2 << endl;
    problems = true;
  }

  v2 = v1 + .3 * wj;

  if ( v1.isNear(v2) ) {
    cout << "OK: with new default tolerance " << v1 <<
	"\n                isNear " << v2 << endl;
  } else {
    cout << "????? with new default tolerance " <<
			v1 << " is not Near " << v2 << endl;
    problems = true;
  }


  v2 = v1 + .45 * wj;

  if ( v1.isNear(v2) ) {
    cout << "????? with new default tolerance " << v1 << " isNear " << v2 << endl;
    problems = true;
  } else {
    cout << "OK: with new default tolerance " << v1 <<
	"\n                is not Near " << v2 << endl;
  }

  v2 = v1 + .6 * wj;

  if ( v1.isNear(v2) ) {
    cout << "?????with new default tolerance " << v1 << " isNear " << v2 << endl;
    problems = true;
  } else {
    cout << "OK: with new default tolerance " << v1 <<
	"\n                is not Near " << v2 << endl;
  }

  v2 = v1 + 1.2 * wj;

  if ( v1.isNear(v2) ) {
    cout << "????? with new default tolerance " << v1 <<
	"\n                isNear " << v2 << endl;
    problems = true;
  } else {
    cout << "OK: " <<
	"with new default tolerance " << v1 <<
	"\n                is not Near " << v2 << endl;
  }

	cout << " { isParallel }\n";

  double percent = .01;
  LorentzVector::setTolerance(percent);
  if ( LorentzVector::getTolerance() == percent ) {
    cout << "OK: new default tolerance is "
		<< LorentzVector::getTolerance() << endl;
  } else {
    cout << "????? new default tolerance is "
		<< LorentzVector::getTolerance() << endl;
    problems = true;
  }

  LorentzVector w1, w2, w3, w4, w5;

  w1.set(2, 4, 5, Tcomponent(8) );
  w2.set(4, 8, 10, Tcomponent(16) );
  w5.set(-10, 4, 8, Tcomponent(18) );

  if ( w1.isParallel(w2,0) ) {
    cout << "OK: " << w1 << " .isParallel " << w2 << " exactly " << endl;
  } else {
    cout << "?????: " << w1 << " .isParallel " << w2 << " FALSE " << endl;
    problems = true;
  }

  w4 = w2;
  w2 = w4 + .3 * percent * w5;
  if ( w1.isParallel(w2) ) {
    cout << "OK: " << w1 << " .isParallel " << w2 << endl;
  } else {
    cout << "?????: " << w1 << " .isParallel\n      " << w4 << " + "
				<< (.3*percent*w5)  << " FALSE " << endl;
    problems = true;
  }

  w2 = w4 + .7 * percent * w5;
  if ( w1.isParallel(w2) ) {
    cout << "OK: " << w1 << " .isParallel " << w2 << endl;
  } else {
    cout << "?????: " << w1 << " .isParallel\n      " << w4 << " + "
				<< (.7*percent*w5)  << " FALSE " << endl;
    problems = true;
  }

  w2 = w4 + 1.1 * percent * w5;
  if ( w1.isParallel(w2) ) {
    cout << "OK: " << w1 << " .isParallel " << w2 << endl;
  } else {
    cout << "?????: " << w1 << " .isParallel\n      " << w4 << " + "
				<< (1.1*percent*w5)  << " FALSE " << endl;
    problems = true;
  }

  // Note that we would expect to find the vectors non-parallel by the time
  // the multiplier used here is about 1.44 times the tolerance.

  w2 = w4 + 1.5 * percent * w5;
  if ( w1.isParallel(w2) ) {
    cout << "?????: " << w1 << " .isParallel\n      " << w4 << " + "
				<< (1.5*percent*w5) << 	endl;
    problems = true;
  } else {
    cout << "OK: NOT " << w1 << " .isParallel " << w2 << endl;
  }

  w2 = w4 + 1.8 * percent * w5;
  if ( w1.isParallel(w2) ) {
    cout << "?????: " << w1 << " .isParallel\n      " << w4 << " + "
				<< (1.8*percent*w5) << 	endl;
    problems = true;
  } else {
    cout << "OK: NOT " << w1 << " .isParallel " << w2 << endl;
  }

} // end of tolerance


void nearCM() {

	cout << " { isNearCM() }\n";

  {
    cout << "  --> part 1:  start near, end up near" << endl;
    double tolerance = 0.000007;

    LorentzVector LV1( SpaceVector( .0001,  .0002,  .0007),  300.    ),
		  LV2( SpaceVector(-.0001, -.0002, -.0007),  299.999 );
    tolerance = 0.000007;
    cout << "Using tolerance of " << tolerance
	 << " with fNear() = " << fNear(LV1, LV2) << ':' << endl;
    if ( LV1.isNear( LV2, tolerance ) )  {
      cout << "OK: " << LV1 << "\n\tisNear " << LV2 << endl;
    } else {
      cout << "????? " << LV1 << "\n\t is not Near " << LV2 << endl;
      problems = true;
    }

    LV1 = boostZOf(LV1, .999999999 );
    LV2 = boostZOf(LV2, .999999999 );
    cout << "After boostZOf(.999999999)"
	 << " with fNear() = " << fNear(LV1, LV2) << ':' << endl;
    if ( ! LV1.isNear( LV2, tolerance ) )  {
      cout << "OK: " << LV1 << "\n\tis not Near " << LV2 << endl;
    } else {
      cout << "????? " << LV1 << "\n\tisNear " << LV2 << endl;
      problems = true;
    }
    if ( LV1.isNearCM( LV2, tolerance ) )  {
      cout << "OK: " << LV1 << "\n\tisNearCM " << LV2 << endl;
    } else {
      cout << "????? " << LV1 << "\n\tis not NearCM " << LV2 << endl;
      problems = true;
    }
  }

  {
    cout << "  --> part 2:  start not near, end up not near" << endl;
    double tolerance = 0.000003;

    LorentzVector LV1( SpaceVector( .0001,  .0002,  .0007),  300.    ),
		  LV2( SpaceVector(-.0001, -.0002, -.0007),  299.999 );
    cout << "Using tolerance of " << tolerance
	 << " with fNear() = " << fNear(LV1, LV2) << ':' << endl;
    if ( ! LV1.isNear( LV2, tolerance ) )  {
      cout << "OK: " << LV1 << "\n\tis not Near " << LV2 << endl;
    } else {
      cout << "????? " << LV1 << "\n\t isNear " << LV2 << endl;
      problems = true;
    }

    LV1 = boostZOf(LV1, -.999999999 );
    LV2 = boostZOf(LV2, -.999999999 );
    cout << "After boostZOf(-.999999999)"
	 << " with fNear() = " << fNear(LV1, LV2) << ':' << endl;
    if ( LV1.isNear( LV2, tolerance ) )  {
      cout << "OK: " << LV1 << "\n\tisNear " << LV2 << endl;
    } else {
      cout << "????? " << LV1 << "\n\tis not Near " << LV2 << endl;
      problems = true;
    }
    if ( ! LV1.isNearCM( LV2, tolerance ) )  {
      cout << "OK: " << LV1 << "\n\tis not NearCM " << LV2 << endl;
    } else {
      cout << "????? " << LV1 << "\n\tisNearCM " << LV2 << endl;
      problems = true;
    }
  }

  {
    cout << "  --> part 3:  validate isNear on vectors in the CM" << endl;
    double tolerance = 0.00003;

    LorentzVector LV1( SpaceVector( .0001,  .0002,  .0007),  100. ),
		  LV2( SpaceVector(-.0001, -.0002, -.0007),  100. );
    cout << "Using tolerance of " << tolerance
	 << " with fNear() = " << fNear(LV1, LV2) << ':' << endl;
    if ( LV1.isNear( LV2, tolerance )  ==  LV1.isNearCM( LV2, tolerance) )  {
      cout << "OK: " << LV1 << "\n\t is"
	   << (LV1.isNear(LV2, tolerance)? "" : " not ") << "Near " << " and is"
	   << (LV1.isNearCM(LV2, tolerance)? "" : " not ") << "NearCM "
	   << LV2 << endl;
    } else {
      cout << "????? " << LV1 << "\n\t is"
	   << (LV1.isNear(LV2, tolerance)? "" : " not ") << "Near " << " but is"
	   << (LV1.isNearCM(LV2, tolerance)? "" : " not ") << "NearCM "
	   << LV2 << endl;
      problems = true;
    }
  }

	cout << " { howNearCM() }\n";

  LorentzVector w6 (  1,  2,  3, Tcomponent (60) );
  LorentzVector w7 ( -1, -2, -3, Tcomponent (64) );

  if (  close (w6.howNearCM(w7),sqrt(72./3858)) &&
	close (w6.howNear(w7),sqrt(72./3858)) ) {
    cout << "OK: " << w6 << "     .howNearCM" << w7 << "= "
		<< w6.howNearCM(w7) << endl;
  } else {
    cout << "?????" << w6 << "     .howNearCM" << w7 << "= "
		<< w6.howNearCM(w7) << endl
		<< w6 << "     .howNear" << w7 << "= "
		<< w6.howNear(w7) << endl;
    problems = true;
  }

  UnitVector bd ( 15, DEGREES, 30, DEGREES );
  w6.boost ( bd, .995 );
  w7.boost ( bd, .995 );
  if (  close (w6.howNearCM(w7),sqrt(72./3858)) ) {
    cout << "OK: " << w6 << "     .howNearCM" << w7 << "= "
		<< w6.howNearCM(w7) << endl;
  } else {
    cout << "?????" << w6 << "     .howNearCM" << w7 << "= "
		<< w6.howNearCM(w7) << endl;
    problems = true;
  }

} // end of nearCM()

void nearness() {

  cout << "\n --- Nearness methods \n\n";

  LorentzVector w1 ( 1, 2, 3, Tcomponent (4) );
  LorentzVector w2 ( 1, 2, 5, Tcomponent (6) );

  if ( close (w1.howNear(w2), sqrt(8./45)) ) {
    cout << "OK: " << w1 << ".howNear" << w2 << "= " << w1.howNear(w2) << endl;
  } else {
    cout << "????? " << w1 << ".howNear" << w2 << "= " << w1.howNear(w2) << endl;
   problems = true;
  }

  w1.setZ( -3 );
  if ( w1.howNear(w2) == 1 ) {
    cout << "OK: " << w1 << ".howNear" << w2 << "= " << w1.howNear(w2) << endl;
  } else {
    cout << "????? " << w1 << ".howNear" << w2 << "= " << w1.howNear(w2) << endl;
   problems = true;
  }

  w1.set (1, 2, -1, Tcomponent(4) );
  w2.set (2, 4,  0, Tcomponent(7) );

  if ( close (w1.howParallel(w2), sqrt(2-76/sqrt(22.*69))  ) ) {
    cout << "OK: " << w1 << ".howParallel" << w2 << "= "
		<< w1.howParallel(w2) << endl;
  } else {
    cout << "????? " << w1 << ".howParallel" << w2 << "= "
		<< w1.howParallel(w2) << endl;
   problems = true;
  }

  w2.set (2, -4,  0, Tcomponent(7) );
  if ( close (w1.howParallel(w2), sqrt(2-44/sqrt(22.*69))  ) ) {
    cout << "OK: " << w1 << ".howParallel" << w2 << "= "
		<< w1.howParallel(w2) << endl;
  } else {
    cout << "????? " << w1 << ".howParallel" << w2 << "= "
		<< w1.howParallel(w2) << endl;
   problems = true;
  }

  w2.set (-4, -2,  0, Tcomponent(7) );
  if ( close (w1.howParallel(w2), sqrt(2-40/sqrt(22.*69))  ) ) {
    cout << "OK: " << w1 << ".howParallel" << w2 << "= "
		<< w1.howParallel(w2) << endl;
  } else {
    cout << "????? " << w1 << ".howParallel" << w2 << "= "
		<< w1.howParallel(w2) << endl;
    cout << "Discrepancy is "
		<< w1.howParallel(w2) - sqrt(2-40/sqrt(22.*69)) << endl;
   problems = true;
  }

  w2.set (-7, -2,  0, Tcomponent(4) );
  if ( w1.howParallel(w2) == 1 ) {
    cout << "OK: " << w1 << ".howParallel" << w2 << "= "
		<< w1.howParallel(w2) << endl;
  } else {
    cout << "????? " << w1 << ".howParallel" << w2 << "= "
		<< w1.howParallel(w2) << endl;
   problems = true;
  }

  LorentzVector w3 ( SpaceVector (10, 1.2, ETA,  .25, RADIANS) , 20 );
  LorentzVector w4 ( SpaceVector (25, 0.9, ETA, -.15, RADIANS) , 30 );

  if ( close(w3.deltaR(w4), .5) ) {
    cout << "OK: " << w3 << "\n    .deltaR" << w4 << "= "
		<< w3.deltaR(v4) << endl;
  } else {
    cout << "????? " << w3 << "\n    .deltaR" << w4 << "= "
		<< w3.deltaR(w4) << endl;
   problems = true;
  }

} // nearness()


void combining() {

// Combining vectors

  cout << "\n --- Methods combining two 4-vectors\n\n";

  if (LorentzVector::getMetric()==TimePositive) {
    cout << "OK: metric is - - - +\n";
  } else {
    cout << "????? metric is + + + -\n";
    problems = true;
  }

	cout << " { dot product }\n";

  v1.set( 3, 5, 1, Tcomponent(30) );
  v2.set( Tcomponent(4), 4, 8, 7);
  double xx;
  xx = v1.dot(v2);
  if (xx==61) {
    cout << "OK: " << v1 << ".dot " << v2 << " = " << xx << endl;
  } else {
    cout << "????? " << v1 << ".dot " << v2 << " = " << xx << endl;
    problems = true;
  }

  v1.set(2, 4, 5, Tcomponent(20));
  v2.set(6, 1, -4, Tcomponent(3));

  double vdv = v1.dot(v2);
  if (vdv == 64) {
    cout << "OK: " << v1 << ".dot " << v2 << " ==> " << vdv << endl;
  } else {
    cout << "????? " << v1 << ".dot " << v2 << " ==> " << vdv << endl;
    problems = true;
  }

  cout << "    ... setting metric to + + + -\n";
  LorentzVector::setMetric(TimeNegative);

  vdv = v2.dot(v1);
  if (vdv == -64) {
    cout << "OK: " << v2 << ".dot " << v1 << " ==> " << vdv << endl;
  } else {
    cout << "????? " << v2 << ".dot " << v1 << " ==> " << vdv << endl;
    problems = true;
  }

  cout << "    ... setting metric back to - - - +\n";
  LorentzVector::setMetric(TimePositive);

  vdv = v2.dot(v2);
  if (vdv == -44) {
    cout << "OK: " << v2 << ".dot " << v2 << " ==> " << vdv << endl;
  } else {
    cout << "????? " << v2 << ".dot " << v2 << " ==> " << vdv << endl;
    problems = true;
  }

	cout << " { Euclidian norm difference squared }\n";

  Float8 vdiff2 = v1.delta2Euclidean(v2);
  if ( vdiff2 ==  395 ) {
    cout << "OK: " << v1 << " delta2Euclidean " << v2 << " ==> "
		<< vdiff2 << endl;
  } else {
    cout << "????? " << v1 << " delta2Euclidean " << v2 << " ==> "
		<< vdiff2 << endl;
    problems = true;
  }

} // end of combining()

void intrinsic() {

  cout << "\n --- Intrinsic properties of a 4-vector\n\n";

  cout << "    ... Setting tolerance to .00005\n";
  LorentzVector::setTolerance (.00005);

  v1.setV( vx );

  v1.setT( 13.001 );

  if ( v1.isTimelike() ) {
    cout << "OK: " << v1 << " isTimelike" << endl;
  } else {
    cout << "????? " << v1 << " is not Timelike" << endl;
    problems = true;
  }

  if ( !v1.isSpacelike() ) {
    cout << "OK: " << v1 << " is not Spacelike" << endl;
  } else {
    cout << "????? " << v1 << " isSpacelike" << endl;
    problems = true;
  }

  if ( !v1.isLightlike() ) {
    cout << "OK: " << v1 << " is not Lightlike" << endl;
  } else {
    cout << "????? " << v1 << " isLightlike" << endl;
    problems = true;
  }

  if ( v1.isLightlike(.0001) ) {
    cout << "OK: " << v1 << " is Lightlike to tolerance .0001" << endl;
  } else {
    cout << "????? " << v1 << " is not Lightlike to tolerance .0001" << endl;
    problems = true;
  }

  double ll = .026001 / ( 338.052002 );
  if ( close(v1.howLightlike(), ll) ) {
    cout << "OK: " << v1 << ".howLightlike() = " << v1.howLightlike() << endl;
  } else {
    cout << "?????" << v1 << ".howLightlike() = " << v1.howLightlike() << endl;
    problems = true;
  }

  v1.setT( 12.999 );

  if ( v1.isSpacelike() ) {
    cout << "OK: " << v1 << " isSpacelike" << endl;
  } else {
    cout << "????? " << v1 << " is not Spacelike" << endl;
    problems = true;
  }

  if ( !v1.isTimelike() ) {
    cout << "OK: " << v1 << " is not Timelike" << endl;
  } else {
    cout << "????? " << v1 << " isTimelike" << endl;
    problems = true;
  }

  if ( !v1.isLightlike() ) {
    cout << "OK: " << v1 << " is not Lightlike" << endl;
  } else {
    cout << "????? " << v1 << " isLightlike" << endl;
    problems = true;
  }

  if ( v1.isLightlike(.0001) ) {
    cout << "OK: " << v1 << " is Lightlike to tolerance .0001" << endl;
  } else {
    cout << "????? " << v1 << " is not Lightlike to tolerance .0001" << endl;
    problems = true;
  }

  cout << "    ... setting metric to + + + -\n";
  LorentzVector::setMetric(TimeNegative);

  v1.setT( 13.001 );

  if ( v1.isTimelike() ) {
    cout << "OK: " << v1 << " isTimelike" << endl;
  } else {
    cout << "????? " << v1 << " is not Timelike" << endl;
    problems = true;
  }

  if ( !v1.isSpacelike() ) {
    cout << "OK: " << v1 << " is not Spacelike" << endl;
  } else {
    cout << "????? " << v1 << " isSpacelike" << endl;
    problems = true;
  }

  if ( !v1.isLightlike() ) {
    cout << "OK: " << v1 << " is not Lightlike" << endl;
  } else {
    cout << "????? " << v1 << " isLightlike" << endl;
    problems = true;
  }

  if ( v1.isLightlike(.0001) ) {
    cout << "OK: " << v1 << " is Lightlike to tolerance .0001" << endl;
  } else {
    cout << "????? " << v1 << " is not Lightlike to tolerance .0001" << endl;
    problems = true;
  }

  v1.setT( 12.999 );

  if ( v1.isSpacelike() ) {
    cout << "OK: " << v1 << " isSpacelike" << endl;
  } else {
    cout << "????? " << v1 << " is not Spacelike" << endl;
    problems = true;
  }

  if ( !v1.isTimelike() ) {
    cout << "OK: " << v1 << " is not Timelike" << endl;
  } else {
    cout << "????? " << v1 << " isTimelike" << endl;
    problems = true;
  }

  if ( !v1.isLightlike() ) {
    cout << "OK: " << v1 << " is not Lightlike" << endl;
  } else {
    cout << "????? " << v1 << " isLightlike" << endl;
    problems = true;
  }

  if ( v1.isLightlike(.0001) ) {
    cout << "OK: " << v1 << " is Lightlike to tolerance .0001" << endl;
  } else {
    cout << "????? " << v1 << " is not Lightlike to tolerance .0001" << endl;
    problems = true;
  }

  v1.setT( 6 );

  if (v1.mag2() == 133) {
    cout << "OK: mag2" << v1 << " = " << v1.mag2() << endl;
  } else {
    cout << "????? mag2" << v1 << " = " << v1.mag2() << endl;
    problems = true;
  }

  cout << "    ... setting metric back to - - - +\n";
  LorentzVector::setMetric(TimePositive);

  if (v1.mag2() == -133) {
    cout << "OK: mag2" << v1 << " = " << v1.mag2() << endl;
  } else {
    cout << "????? mag2" << v1 << " = " << v1.mag2() << endl;
    problems = true;
  }

  if (v1.plus() == 18) {
    cout << "OK: " << v1 << ".plus() = " << v1.plus() << endl;
  } else {
    cout << "????? " << v1 << ".plus() = " << v1.plus() << endl;
    problems = true;
  }

  if (v1.minus() == -6) {
    cout << "OK: " << v1 << ".minus() = " << v1.minus() << endl;
  } else {
    cout << "????? " << v1 << ".minus() = " << v1.minus() << endl;
    problems = true;
  }

  UnitVector ux ( 0, -4, -3 );

  if ( close(v1.plus(ux), 2) ) {
    cout << "OK: " << v1 << ".plus " << ux << " = " << v1.plus(ux) << endl;
  } else {
    cout << "????? " << v1 << ".plus " << ux << " = " << v1.plus(ux) << endl;
    problems = true;
  }

  if ( close(v1.minus(ux), 10) ) {
    cout << "OK: " << v1 << ".minus " << ux << " = " << v1.minus(ux) << endl;
  } else {
    cout << "????? " << v1 << ".minus " << ux << " = " << v1.minus(ux) << endl;
    problems = true;
  }

  vx.set ( 0, -4, -3 );

  if ( close(v1.plus(vx), 2) ) {
    cout << "OK: " << v1 << ".plus " << vx << " = " << v1.plus(vx) << endl;
  } else {
    cout << "????? " << v1 << ".plus " << vx << " = " << v1.plus(vx) << endl;
    problems = true;
  }

  if ( close(v1.minus(vx), 10) ) {
    cout << "OK: " << v1 << ".minus " << vx << " = " << v1.minus(vx) << endl;
  } else {
    cout << "????? " << v1 << ".minus " << vx << " = " << v1.minus(vx) << endl;
    problems = true;
  }

  v1.set ( 48, -12, 16, Tcomponent(39) );

  if ( v1.euclideanNorm2() == 65*65 ) {
    cout << "OK: " << v1 << ".euclideanNorm2() = "
				<<  v1.euclideanNorm2() << endl;
  } else {
    cout << "????? " << v1 << ".euclideanNorm2() = "
				<<  v1.euclideanNorm2() << endl;
    problems = true;
  }

  if ( v1.euclideanNorm() == 65) {
    cout << "OK: " << v1 << ".euclideanNorm() = "
				<<  v1.euclideanNorm() << endl;
  } else {
    cout << "????? " << v1 << ".euclideanNorm() = "
				<<  v1.euclideanNorm() << endl;
    problems = true;
  }

} // end of intrinsic()

void kinematic() {

// Relativistic kinematic properties

  cout << "\n --- Relativistic kinematic properties\n\n";

	cout << " { rest-frame quantities }\n";

  v1.set ( Tcomponent(65), 48, 12, -16);

  if ( v1.restMass2() == 39*39 ) {
    cout << "OK: " << v1 << ".restMass2() = "
				<<  v1.restMass2() << endl;
  } else {
    cout << "????? " << v1 << ".restMass2() = "
				<<  v1.restMass2() << endl;
    problems = true;
  }

// THE FOLLOWING JUST TESTS EXCEPTIONS AND SHOULD BE REMOVED
//  v1.set ( Tcomponent(1), 48, 12, -16);
//  if ( v1.restMass2() == 39*39 ) {
//    cout << "OK: " << v1 << ".restMass2() = "
//				<<  v1.restMass2() << endl;
//  } 

  if ( v1.restMass() == 39 ) {
    cout << "OK: " << v1 << ".restMass() = "
				<<  v1.restMass() << endl;
  } else {
    cout << "????? " << v1 << ".restMass() = "
				<<  v1.restMass() << endl;
    problems = true;
  }

  v1.setT( v1.t() * -1 );

  if ( v1.restMass() == -39 ) {
    cout << "OK: " << v1 << ".restMass() = "
				<<  v1.restMass() << endl;
  } else {
    cout << "????? " << v1 << ".restMass() = "
				<<  v1.restMass() << endl;
    problems = true;
  }

  v1.set ( Tcomponent(65), 48, 12, -16);

  if (eq(v1.beta() , 4./5 ) ) {
    cout << "OK: " << v1 << ".beta() == " << v1.beta() << endl;
  } else {
    cout << "????? " << v1 << ".beta() == " << v1.beta() << endl;
    problems = true;
  }

  if ( close(v1.gamma(), 5./3) ) {
    cout << "OK: " << v1 << ".gamma() == " << v1.gamma() << endl;
  } else {
    cout << "????? " << v1 << ".gamma() == " << v1.gamma() << endl;
    problems = true;
  }

  LorentzVector vm ( Tcomponent (39) );
  LorentzVector vmcopy ( vm );

  if ( v1.rest4Vector() == vm ) {
    cout << "OK: " << v1 << ".rest4Vector() = " << v1.rest4Vector() << endl;
  } else {
    cout << "????? " << v1 << ".rest4Vector() = " << v1.rest4Vector() << endl;
    problems = true;
  }

  cout << "    ... Non-trivial test of boost from rest frame:\n";

  LorentzVector::setTolerance(toler);

  vx = v1.boostVector();
  vm.boost(vx);
  if (vm.isNear(v1)) {
    cout << "OK: " << vmcopy << ".boost" << vx <<
	"\n           =" << vm << endl;
  } else {
    cout << "????? " << vmcopy << ".boost" << vx << "=" << vm << endl;
    problems = true;
  }

  v1.set ( Tcomponent(-65), 48, 12, -16);

  if (eq(v1.beta() , 4./5 ) ) {
    cout << "OK: " << v1 << ".beta() == " << v1.beta() << endl;
  } else {
    cout << "????? " << v1 << ".beta() == " << v1.beta() << endl;
    problems = true;
  }

  if ( close(v1.gamma(), 5./3) ) {
    cout << "OK: " << v1 << ".gamma() == " << v1.gamma() << endl;
  } else {
    cout << "????? " << v1 << ".gamma() == " << v1.gamma() << endl;
    problems = true;
  }

  vm.set ( (-39) );
  vmcopy = vm;

  if ( v1.rest4Vector() == vm ) {
    cout << "OK: " << v1 << ".rest4Vector() = " << v1.rest4Vector() << endl;
  } else {
    cout << "????? " << v1 << ".rest4Vector() = " << v1.rest4Vector() << endl;
    problems = true;
  }

  cout << "    ... Non-trivial test of boost from rest frame:\n";
  vx = v1.boostVector();
  vm.boost(vx);
  if (vm.isNear(v1)) {
    cout << "OK: " << vmcopy << ".boost" << vx <<
	"\n                 =" << vm << endl;
  } else {
    cout << "????? " << vmcopy << ".boost" << vx <<
	"\n                =" << vm << endl;
    problems = true;
  }


  cout << "    ... Transverse energy:\n";

  v1.set ( Tcomponent(52.000000001), -12, 16, -48 );  // small mass (.00032)
  double Et = v1.et();
  cout << "      " << v1 << " et() = " << Et << "\n";
  vx.set ( 0, 0, .6 );
  double EtPrime = (boostOf(v1,vx)).et();  
  // Et^2 - EtPrime^2 should be pperp^2 m^2 (1/E^2 - 1/Eprime^2)
  double Eprime = (boostOf(v1,vx)).e();
  double Eprime2 = Eprime*Eprime;
  double E2 = v1.e()*v1.e();
  double first_order = v1.perp2() * v1.m2() * ( 1/E2 - 1/Eprime2 );
  double difference = Et*Et - EtPrime*EtPrime;
  cout << "      " << (boostOf(v1,vx)) << " et() = " << EtPrime << "\n";
  cout << "      " << "difference of squares is " << difference << 
	  "; estimate is " << first_order << "\n";
  if ( close ( difference, first_order ) ) {
    cout << "OK: " << v1 << "transverse energy transforms as expected\n";
  } else {
    cout << "????? " << v1 << "transverse energy transforms improperly\n";
    problems = true;
  }

  vx.set ( .3, .4, .1 );
  Et = v1.et(vx);
  cout << "      "<< v1 << " et(vx) = " << Et << "\n";
  EtPrime = (boostOf(v1,.6*vx.unit())).et(vx);  
  // Et^2 - EtPrime^2 should be pperp^2 m^2 (1/E^2 - 1/Eprime^2)
  Eprime = (boostOf(v1,.6*vx.unit())).e();
  Eprime2 = Eprime*Eprime;
  E2 = v1.e()*v1.e();
  first_order = v1.perp2(vx) * v1.m2() * ( 1/E2 - 1/Eprime2 );
  difference = Et*Et - EtPrime*EtPrime;
  cout << "      "<< (boostOf(v1,.6*vx.unit())) << ".et" << vx << " = " 
			<< EtPrime << "\n";
  cout << "      "<< "difference of squares is " << difference << 
	  "; estimate is " << first_order << "\n";
  if ( close ( difference, first_order ) ) {
    cout << "OK: " << v1 << ".et" << vx << "transforms as expected\n";
  } else {
    cout << "????? " << v1 << ".et" << vx << "transforms improperly\n";
    problems = true;
  }


	cout << " { pseudorapidity and rapidity }\n";

  v1.set ( Tcomponent(65), -48, -12, 16);
  if ( eq(v1.eta(), v1.getV().eta()) ) {
    cout << "OK: " << v1 << ".eta() = " << v1.eta() << endl;
  } else {
    cout << "????? " << v1 << ".eta() = " << v1.eta() << endl;
    problems = true;
  }

  vx.set( 1, 2, 3);
  if ( close(v1.eta(vx), v1.getV().eta(UnitVector(vx))) ) {
    cout << "OK: " << v1 << ".eta() w.r.t. " << UnitVector(vx) <<
	"\n      = " << v1.eta() << endl;
  } else {
    cout << "?????" << v1 << ".eta() w.r.t. " << UnitVector(vx) <<
	"\n      = " << v1.eta() << endl;
    problems = true;
  }

  cout << "    ... Non-trivial check that tanh(rapidity) = z/t\n";
  Float8 rap = v1.rapidity();
  Float8 tanhrap =  (exp(rap) - exp(-rap)) / (exp(rap) + exp(-rap));
  if ( close(tanhrap, v1.z()/v1.t()) ) {
    cout << "OK:" << v1 << ".rapidity() = " << rap
	<< "\n              tanh(rapidity) = " << tanhrap
	<< " z/t = " << v1.z()/v1.t() << endl;
  } else {
    cout << "?????" << v1 << ".rapidity() = " << rap
	<< "\n              tanh(rapidity) = " << tanhrap
	<< " z/t = " << v1.z()/v1.t() << endl;
    problems = true;
  }

  cout << "    ... Non-trivial check that tanh(w.rapidity(u)) = v.dot(u)/t\n";
  UnitVector ux = UnitVector(vx);
  rap = v1.rapidity(ux);
  tanhrap =  (exp(rap) - exp(-rap)) / (exp(rap) + exp(-rap));
  if ( close(tanhrap, v1.getV().dot(ux)/v1.t()) ) {
    cout << "OK:" << v1 << ".rapidity(u) = " << rap
	<< "\n              tanh(rapidity(u)) = " << tanhrap
	<< " v.dot(u)/t = " << v1.getV().dot(ux)/v1.t() << endl;
  } else {
    cout << "?????" << v1 << ".rapidity(u) = " << rap
	<< "\n              tanh(rapidity(u)) = " << tanhrap
	<< " v.dot(u)/t = " << v1.getV().dot(ux)/v1.t() << endl;
    problems = true;
  }

  cout << "    ... Non-trivial check that tanh(coLinearRapidity) = v/t\n";
  rap = v1.coLinearRapidity();
  tanhrap =  (exp(rap) - exp(-rap)) / (exp(rap) + exp(-rap));
  if ( close(tanhrap, v1.getV().mag()/v1.t()) ) {
    cout << "OK:" << v1 << ".coLinearRapidity() = " << rap
	<< "\n              tanh(rapidity) = " << tanhrap
	<< " v/t = " << v1.getV().mag()/v1.t() << endl;
  } else {
    cout << "?????" << v1 << ".coLinearRapidity() = " << rap
	<< "\n              tanh(rapidity) = " << tanhrap
	<< " v/t = " << v1.getV().mag()/v1.t() << endl;
    problems = true;
  }

	cout << " { Invariant mass of a pair of particles }\n";

  v1.set ( 6,  7, -25, Tcomponent(29) );
  v2.set ( 10, 5, -23, Tcomponent(36) );

  if ( close (v1.invariantMass2(v2), 39*39) ) {
    cout << "OK: " << v1 << ".invariantMass2" << v2 << " = " <<
		v1.invariantMass2(v2) << endl;
  } else {
    cout << "?????" << v1 << ".invariantMass2" << v2 << " = " <<
		v1.invariantMass2(v2) << endl;
    problems = true;
  }

  if ( close (v2.invariantMass2(v1), 39*39) ) {
    cout << "OK: " << v2 << ".invariantMass2" << v1 << " = " <<
		v2.invariantMass2(v1) << endl;
  } else {
    cout << "?????" << v2 << ".invariantMass2" << v1 << " = " <<
		v2.invariantMass2(v1) << endl;
    problems = true;
  }

  if ( close (v1.invariantMass(v2), 39) ) {
    cout << "OK: " << v1 << ".invariantMass " << v2 << " = " <<
		v1.invariantMass(v2) << endl;
  } else {
    cout << "?????" << v1 << ".invariantMass " << v2 << " = " <<
		v1.invariantMass(v2) << endl;
    problems = true;
  }

  if ( close (v2.invariantMass(v1), 39) ) {
    cout << "OK: " << v2 << ".invariantMass " << v1 << " = " <<
		v2.invariantMass(v1) << endl;
  } else {
    cout << "?????" << v2 << ".invariantMass " << v1 << " = " <<
		v2.invariantMass(v1) << endl;
    problems = true;
  }

  if ( close (v1.invariantMass(v1),2*v1.restMass()) ) {
    cout << "OK: " << v1 << ".invariantMass(self) = " << v1.invariantMass(v1)
	<< "\n         twice v1.restMass = " << 2*v1.restMass() << endl;
  } else {
    cout << "????? " << v1 << ".invariantMass(self) = " << v1.invariantMass(v1)
	<< "\n        twice v1.restMass = " << 2*v1.restMass() << endl;
    problems = true;
  }

  cout << "    ... Non-trivial validation that invariantMass is "
	  "invariant after boost\n";

  vx.set (.3, -.4, .2);
  v1.boost(vx);
  v2.boost(vx);
  cout << "    ... Boosting previous vectors by " << vx << endl;

  if ( close (v2.invariantMass(v1), 39) ) {
    cout << "OK: " << v2 << ".invariantMass\n    " << v1 << " = " <<
		v2.invariantMass(v1) << endl;
  } else {
    cout << "?????" << v2 << ".invariantMass\n    " << v1 << " = " <<
		v2.invariantMass(v1) << endl;
    problems = true;
  }

	cout << " { Boost needed to go into CM frame }\n";

  v1.set ( Tcomponent(65), -48, -12, 16);
  vx = v1.findBoostToCM();
  if ( vx.isNear(SpaceVector ( 48./65, 12./65, -16./65 ) ) ) {
    cout << "OK: " << v1 << ".findBoostToCM()" <<
	"\n             =" << vx << "= 1/65" << 65*vx << endl;
  } else {
    cout << "?????" << v1<< ".findBoostToCM()" <<
	"\n             =" << vx << "!= "<<
	SpaceVector ( 48./65, 12./65, -16./65 ) << endl;
    problems = true;
  }

  cout << "    ... Checking single vector Boost-to-CM identities\n";
  if ( v1.findBoostToCM() == - v1.boostVector() ) {
    cout << "OK: " << v1 << ".findBoostToCM() = - v1.boostVector()\n";
  } else {
    cout << "?????" << v1 << ".findBoostToCM() = " << v1.findBoostToCM() <<
	"\n          - v1.boostVector() =" << - v1.boostVector() << endl;
    problems = true;
  }

  v4 = v3 = v2 = v1;
  v2.boost(v2.findBoostToCM());
  v4.boost(v1.findBoostToCM());
  if ( v2.isNear(v1.rest4Vector()) ) {
    cout << "OK: " << v1 << ".boost(v1.findBoostToCM())" <<
	"\n    ="  << v4 << "rest4Vector is" <<
	v3.rest4Vector() << endl;
  } else {
    cout << "?????" << v1 << ".boost(v1.findBoostToCM())" <<
	"\n    ="  << v4 << "rest4Vector is" <<
	v3.rest4Vector() << endl;
    problems = true;
  }

  v1.set ( Tcomponent(15), -14, 0, 9);
  v2.set ( Tcomponent(50), 2, -48, 7);

  vx = v1.findBoostToCM(v2);
  if ( vx.isNear( SpaceVector ( 12./65, 48./65, -16./65 ) ) ) {
    cout << "OK: " << v1 << ".findBoostToCM" << v2 <<
	"\n             =" << vx << "= 1/65" << 65*vx << endl;
  } else {
    cout << "OK: " << v1<< ".findBoostToCM" << v2 <<
	"\n             =" << vx << "!= "<<
	SpaceVector ( 12./65, 48./65, -16./65 ) << endl;
    problems = true;
  }
  cout << "    ... Validating that this really boosts to CM frame:\n";
  if ( boostOf(v1,vx).getV().isNear(-boostOf(v2,vx).getV(),.000000001) ) {
    cout << "OK: " << v1 << "-->" << boostOf(v1,vx) << endl <<
            "    " << v2 << "-->" << boostOf(v2,vx) << endl;
  } else {
    cout << "???? " << v1 << "-->" << boostOf(v1,vx) << endl <<
            "     " << v2 << "-->" << boostOf(v2,vx) << endl;
    problems = true;
  }

} // end of kinematic

void rotations() {

// ROTATIONS

  cout << "\n --- Rotations\n\n";

	cout << " { 90 degree rotations about X Y or Z axis }\n";

  v1.set ( 9, 5, -3, Tcomponent(1) );

  LorentzVector::setTolerance ( .0001 );

  v2 = v1;
  v2.rotateX( PI/2 );
  v3.set ( 9, 3, 5, Tcomponent(1) );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v1 << "rotateX(90 degrees) ="
		"\n                      "   << v2 << endl;
  } else {
    cout << "????? " << v1 << "rotateX(90 degrees) ="
		"\n                      "   << v2 << endl;
    problems=true;
  }
  v2 = v3;
  v2.rotateX( -PI/2 );
  if ( v2.isNear(v1) ) {
    cout << "OK: " << v3 << "rotateX( -PI/2 )    ="
		"\n                      "   << v2 << endl;
  } else {
    cout << "????? " << v3 << "rotateX( -PI/2 ) ="
		"\n                      "   << v2 << endl;
    problems=true;
  }

  v2 = v1;
  v2.rotateY( PI/2 );
  v3.set ( -3, 5, -9, Tcomponent(1)  );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v1 << "rotateY(90 degrees) ="
		"\n                      "   << v2 << endl;
  } else {
    cout << "????? " << v1 << "rotateY(90 degrees) ="
		"\n                      "   << v2 << endl;
    problems=true;
  }
  v2 = v3;
  v2.rotateY( -PI/2 );
  if ( v2.isNear(v1) ) {
    cout << "OK: " << v3 << "rotateY( -PI/2 )    ="
		"\n                      "   << v2 << endl;
  } else {
    cout << "????? " << v3 << "rotateY( -PI/2 ) ="
		"\n                      "   << v2 << endl;
    problems=true;
  }

  v2 = v1;
  v2.rotateZ( PI/2 );
  v3.set ( -5, 9, -3, Tcomponent(1)  );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v1 << "rotateZ(90 degrees) ="
		"\n                      "   << v2 << endl;
  } else {
    cout << "????? " << v1 << "rotateZ(90 degrees) ="
		"\n                      "   << v2 << endl;
    problems=true;
  }
  v2 = v3;
  v2.rotateZ( -PI/2 );
  if ( v2.isNear(v1) ) {
    cout << "OK: " << v3 << "rotateZ( -PI/2 )    ="
		"\n                      "   << v2 << endl;
  } else {
    cout << "????? " << v3 << "rotateZ( -PI/2 ) ="
		"\n                      "   << v2 << endl;
    problems=true;
  }


	cout << " { Pair of 90 degree rotations to permute x y z }\n";

  v1.set ( 9, 5, -3, Tcomponent(1)  );
  v3.set ( -3, 9, 5, Tcomponent(1)  );
  v2 = v1;
  v2.rotateX ( PI/2 );
  v2.rotateZ ( PI/2 );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v1 << "becomes " << v2 << endl;
  } else {
    cout << "????? " << v1 << "becomes " << v2 << endl;
    problems=true;
  }

  v3 = v2;
  v2.rotateZ ( -PI/2 );
  v2.rotateX ( -PI/2 );
  if ( v2.isNear(v1) ) {
    cout << "OK: " << v3 << "becomes " << v2 << endl;
  } else {
    cout << "????? " << v3 << "becomes " << v2 << endl;
    problems=true;
  }


	cout << " { Axis rotations by general angles }\n";

  v1.set ( 1, 10*sqrt(3.), 10, Tcomponent(1)  );
  v3.set ( 1, 0, 20, Tcomponent(1)  );
  v2 = v1;
  v2.rotateX ( PI/3 );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v1 << "rotateX(60 degrees) ="
		"\n                      "   << v2 << endl;
  } else {
    cout << "????? " << v1 << "rotateX(60 degrees) ="
		"\n                      "   << v2 << endl;
    problems=true;
  }
  v2.rotateX ( 2*PI/3 );
  v1.setX( v1.x() * -1 );
  v1.setT( v1.t() * -1 );
  if ( v2.isNear(-v1) ) {
    cout << "OK: rotated another 120 degrees gets back to " << v2 << endl;
  } else {
    cout << "????? rotate another 120 degrees gets back to "
								<< v2 << endl;
    problems=true;
  }

  v1.set ( 10, 1, 10*sqrt(3.), Tcomponent(1)  );
  v3.set ( 20, 1, 0, Tcomponent(1)  );
  v2 = v1;
  v2.rotateY ( PI/3 );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v1 << "rotateY(60 degrees) ="
		"\n                      "   << v2 << endl;
  } else {
    cout << "????? " << v1 << "rotateY(60 degrees) ="
		"\n                      "   << v2 << endl;
    problems=true;
  }
  v2.rotateY ( 2*PI/3 );
  v1.setY( v1.y() * -1 );
  v1.setT( v1.t() * -1 );
  if ( v2.isNear(-v1) ) {
    cout << "OK: rotated another 120 degrees gets back to " << v2 << endl;
  } else {
    cout << "????? rotate another 120 degrees gets back to "
								<< v2 << endl;
    problems=true;
  }

  v1.set ( 10*sqrt(3.), 10, 1, Tcomponent(1)  );
  v3.set ( 0, 20, 1, Tcomponent(1)  );
  v2 = v1;
  v2.rotateZ ( PI/3 );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v1 << "rotateZ(60 degrees) ="
		"\n                      "   << v2 << endl;
  } else {
    cout << "????? " << v1 << "rotateZ(60 degrees) ="
		"\n                      "   << v2 << endl;
    problems=true;
  }
  v2.rotateZ ( 2*PI/3 );
  v1.setZ( v1.z() * -1 );
  v1.setT( v1.t() * -1 );
  if ( v2.isNear(-v1) ) {
    cout << "OK: rotated another 120 degrees gets back to " << v2 << endl;
  } else {
    cout << "????? rotate another 120 degrees gets back to "
								<< v2 << endl;
    problems=true;
  }

} // end of rotations()

void rotations2() {

	cout << " { Generic rotations using XYZ axes }\n";

  v1.set (8, -5, 3, Tcomponent(1) );

  SpaceVector xv (4,0,0);
  UnitVector  xu (1,0,0);
  v3 = v2 = v1;
  v2.rotate ( xu, 50*PI/180 );
  v3.rotateX ( 50*PI/180 );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v2 << "~" << v3 << endl;
  } else {
    cout << "????? " << v2 << " not near " << v3 << endl;
    problems=true;
  }

  v3 = v2;
  v2.rotate ( xv, -.4 );
  v3.rotateX( -.4 );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v2 << "~" << v3 << endl;
  } else {
    cout << "????? " << v2 << " not near " << v3 << endl;
    problems=true;
  }

  v3 = v2;
  v2.rotate ( -xv, .9 );
  v3.rotateX( -.9 );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v2 << "~" << v3 << endl;
  } else {
    cout << "????? " << v2 << " not near " << v3 << endl;
    problems=true;
  }

  SpaceVector yv (0,5,0);
  UnitVector  yu (0,-1,0);
  v3 = v2 = v1;
  v2.rotate ( yu, 50*PI/180 );
  v3.rotateY ( -50*PI/180 );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v2 << "~" << v3 << endl;
  } else {
    cout << "????? " << v2 << " not near " << v3 << endl;
    problems=true;
  }

  v3 = v2;
  v2.rotate ( yv, -.4 );
  v3.rotateY( -.4 );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v2 << "~" << v3 << endl;
  } else {
    cout << "????? " << v2 << " not near " << v3 << endl;
    problems=true;
  }

  v3 = v2;
  v2.rotate ( -yv, .9 );
  v3.rotateY( -.9 );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v2 << "~" << v3 << endl;
  } else {
    cout << "????? " << v2 << " not near " << v3 << endl;
    problems=true;
  }

  SpaceVector zv (0,0,6);
  UnitVector  zu (0,0,1);
  v3 = v2 = v1;
  v2.rotate ( zu, 50*PI/180 );
  v3.rotateZ ( 50*PI/180 );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v2 << "~" << v3 << endl;
  } else {
    cout << "????? " << v2 << " not near " << v3 << endl;
    problems=true;
  }

  v3 = v2;
  v2.rotate ( zv, -.4 );
  v3.rotateZ( -.4 );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v2 << "~" << v3 << endl;
  } else {
    cout << "????? " << v2 << " not near " << v3 << endl;
    problems=true;
  }

  v3 = v2;
  v2.rotate ( -zv, .9 );
  v3.rotateZ( -.9 );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v2 << "~" << v3 << endl;
  } else {
    cout << "????? " << v2 << " not near " << v3 << endl;
    problems=true;
  }


	cout << " { Generic rotations about arbitrary axes }\n";

  // Laboriously hand-calculated :

  UnitVector vu ( 1, 1, 0);

  v1.set ( 3, 4, 0, Tcomponent(1)  );
  v3 = v2 = v1;
  v2.rotate ( vu, PI/4 );
  cout << "    ... " << v1 << "rotated 45 degrees about "
	"\n        " << vu << " is" << 	v2 << endl;
  v3.rotate ( vu, -3*PI/4 );
  cout << "    "<< v1 << "rotated -135 degrees about "
	"\n        " << vu << " is" << 	v3 << endl;
  double temp2 = v3.x(); v3.setX( v3.y() ); v3.setY( temp2 );
  v3.setZ( v3.z() * -1 );
  if ( v2.isNear(v3) ) {
    cout << "OK: " << v2 << "~" << v3 << endl;
  } else {
    cout << "????? " << v2 << " not near " << v3 << endl;
    problems=true;
  }

  // Using formula described below:

  SpaceVector v2s;

  double phi;

  v1.set(3,4,1, Tcomponent(1) );
  v2s.set(1,1,0);
  phi = PI/4;
  v3.set( 3.6464, 3.3536, 1.2071, Tcomponent(1)  );
  if ( checkAxialRotation (v1, v2s, phi, v3 ) ) problems = true;

  v1.set(3,4,1, Tcomponent(1) );
  v2s.set(2,2,0);
  phi = PI/2;
  v3.set( 4.2071, 2.7929, 0.7071, Tcomponent(1)  );
  if ( checkAxialRotation (v1, v2s, phi, v3 ) ) problems = true;

  v1.set(3,4,1, Tcomponent(1) );
  v2s.set(1,1,0);
  phi = 175*PI/180;
  v3.set( 4.0597, 2.9403, -0.9346, Tcomponent(1)  );
  if ( checkAxialRotation (v1, v2s, phi, v3 ) ) problems = true;

#ifdef INTENTIONAL_ERROR
  // Intentionally wrong to check our checker:
  v1.set(3,4,1, Tcomponent(1) );
  v2s.set(1,1,0);
  phi = 175*PI/180;
  v3.set( 4.0397, 2.9403, -0.9346, Tcomponent(1)  );
  if ( checkAxialRotation (v1, v2s, phi, v3 ) ) problems = true;
#endif // INTENTIONAL_ERROR

  v1.set(2,-1,3, Tcomponent(1) );
  v2s.set(4,2,-5);
  phi = 35*PI/180;
  v3.set( 1.5791, -2.7726, 1.9543, Tcomponent(1)  );
  if ( checkAxialRotation (v1, v2s, phi, v3 ) ) problems = true;

  v1.set(2,-1,3, Tcomponent(1) );
  v2s.set(4,2,-5);
  phi = -10*PI/180;
  v3.set( 1.9316, -0.4214, 3.1767, Tcomponent(1)  );
  if ( checkAxialRotation (v1, v2s, phi, v3 ) ) problems = true;

  v1.set(2,-1,3, Tcomponent(1) );
  v2s.set(4,2,-5);
  phi = -100*PI/180;
  v3.set( -1.4330, 2.9340, 1.8272, Tcomponent(1)  );
  if ( checkAxialRotation (v1, v2s, phi, v3 ) ) problems = true;

  v1.set(2,3,4, Tcomponent(1) );
  v2s.set(0,1,0);
  phi = 20*PI/180;
  v3.set( 3.2475, 3, 3.0747, Tcomponent(1) );
  if ( checkAxialRotation (v1, v2s, phi, v3 ) ) problems = true;

  v1.set(2,3,4, Tcomponent(1) );
  v2s.set(0,3,0);
  phi = PI/2;
  v3.set(4, 3, -2, Tcomponent(1) );
  if ( checkAxialRotation (v1, v2s, phi, v3 ) ) problems = true;

  v1.set(2,3,4, Tcomponent(1) );
  v2s.set(0,1,0);
  phi = PI;
  v3.set( -2, 3, -4, Tcomponent(1)  );
  if ( checkAxialRotation (v1, v2s, phi, v3 ) ) problems = true;

  v1.set(2,3,4, Tcomponent(1) );
  v2s.set(0,-1,0);
  phi = 20*PI/180;
  v3.set( 0.5113, 3, 4.4428, Tcomponent(1)  );
  if ( checkAxialRotation (v1, v2s, phi, v3 ) ) problems = true;

  v1.set(2,3,4, Tcomponent(1) );
  v2s.set(0,-1,0);
  phi = -20*PI/180;
  v3.set( 3.2475, 3, 3.0747, Tcomponent(1) );
  if ( checkAxialRotation (v1, v2s, phi, v3 ) ) problems = true;

  v1.set(1,3,4, Tcomponent(1) );
  v2s.set(0,.01,1);
  phi = 1*PI/180;
  v3.set( 0.9482, 3.0170, 3.9998, Tcomponent(1)  );
  if ( checkAxialRotation (v1, v2s, phi, v3 ) ) problems = true;

  v1.set(1,3,4, Tcomponent(1) );
  v2s.set(0,.01,1);
  phi = -10*PI/180;
  v3.set( 1.4988, 2.7814, 4.0022, Tcomponent(1)  );
  if ( checkAxialRotation (v1, v2s, phi, v3 ) ) problems = true;

  v1.set(2,3,4, Tcomponent(1) );
  v2s.set(1, -1, 1);
  phi = -130*PI/180;
  v3.set( 3.4531, -2.6866, -3.1397 , Tcomponent(1) );
  if ( checkAxialRotation (v1, v2s, phi, v3 ) ) problems = true;

  v1.set(2,3,4, Tcomponent(1) );
  v2s.set(1, -1, 1);
  phi = -80*PI/180;
  v3.set( 5.1537, 0.8318, -1.3220, Tcomponent(1)  );
  if ( checkAxialRotation (v1, v2s, phi, v3 ) ) problems = true;

  v1.set(2,3,4, Tcomponent(1) );
  v2s.set(1, -1, 1);
  phi = -PI/6;
  v3.set( 3.8868, 3.0415, 2.1547, Tcomponent(1)  );
  if ( checkAxialRotation (v1, v2s, phi, v3 ) ) problems = true;

  v1.set(2,3,4, Tcomponent(1) );
  v2s.set(1, -1, 1);
  phi = 20*PI/180;
  v3.set( 0.5574, 2.3638, 4.8064, Tcomponent(1)  );
  if ( checkAxialRotation (v1, v2s, phi, v3 ) ) problems = true;

  v1.set(2,3,4, Tcomponent(1) );
  v2s.set(1, -1, 1);
  phi = 70*PI/180;
  v3.set( -2.4557, -0.7170, 4.7387, Tcomponent(1)  );
  if ( checkAxialRotation (v1, v2s, phi, v3 ) ) problems = true;

} // end of rotations2()


void rotations3() {

}  // end of rotations2()


void boosts() {

  cout << "\n --- Boosts\n\n";

  cout << "{ Non-trivial tests of rapidity addition under boosts }\n";

  v1.set   ( Tcomponent(70), 31, -25, 46);
  v2 = boostZOf(v1,.24);
  vx.set(0,0,.24);
  if ( close(v2.rapidity(),v1.rapidity()+vx.coLinearRapidity()) ) {
    cout << "OK: " << v1 << "(rapidity " << v1.rapidity() << ")\n"
        << "    boosted by " << vx
	<< "(rapidity " << vx.coLinearRapidity() << ")\n"
	<< "    =" << v2 << "(rapidity " << v2.rapidity() << " ~ "
	<< v1.rapidity() + vx.coLinearRapidity() << endl;
  } else {
    cout << "?????" << v1 << "(rapidity " << v1.rapidity() << ")\n"
        << "    boosted by " << vx
	<< "(rapidity " << vx.coLinearRapidity() << ")\n"
	<< "    =" << v2 << "(rapidity " << v2.rapidity() << " ~ "
	<< v1.rapidity() + vx.coLinearRapidity() << endl;
    problems = true;
  }

  v1.set   ( Tcomponent(90), 51, 25, 36);
  v2 = boostYOf(v1,.84);
  vx.set(0,.84,0);
  if ( close(v2.rapidity(Y_HAT),v1.rapidity(Y_HAT)+vx.coLinearRapidity()) ) {
    cout << "OK: " << v1 << "(Y-rapidity " << v1.rapidity(Y_HAT) << ")\n"
        << "    boosted by " << vx
	<< "(rapidity " << vx.coLinearRapidity() << ")\n"
	<< "    =" << v2 << "(rapidity " << v2.rapidity(Y_HAT) << " ~ "
	<< v1.rapidity(Y_HAT) + vx.coLinearRapidity() << endl;
  } else {
    cout << "?????" << v1 << "(Y-rapidity " << v1.rapidity(Y_HAT) << ")\n"
        << "    boosted by " << vx
	<< "(rapidity " << vx.coLinearRapidity() << ")\n"
	<< "    =" << v2 << "(rapidity " << v2.rapidity(Y_HAT) << " ~ "
	<< v1.rapidity(Y_HAT) + vx.coLinearRapidity() << endl;
    problems = true;
  }

  cout << "    ... Here rapidity reference direction is supplied "
			"as a SpaceVector:\n";
  v1.set   ( Tcomponent(80), -31, 35, 26);
  v2 = boostXOf(v1,-.53);
  vx.set(-.53,0,0);
		// Note negative X boost:  Since colinear rapidity
		// XXXXXXXXXXXX
  SpaceVector longvx(3,0,0);
  if ( close(v2.rapidity(longvx),v1.rapidity(X_HAT)-vx.coLinearRapidity()) ) {
    cout << "OK: " << v1 << "(X-rapidity " << v1.rapidity(X_HAT) << ")\n"
        << "    boosted by " << vx
	<< "(rapidity " << vx.coLinearRapidity() << ")\n"
	<< "    =" << v2 << "(rapidity " << v2.rapidity(longvx) << " ~ "
	<< v1.rapidity(X_HAT) - vx.coLinearRapidity() << endl;
  } else {
    cout << "?????" << v1 << "(X-rapidity " << v1.rapidity(X_HAT) << ")\n"
        << "    boosted by " << vx
	<< "(rapidity " << vx.coLinearRapidity() << ")\n"
	<< "    =" << v2 << "(rapidity " << v2.rapidity(longvx) << " ~ "
	<< v1.rapidity(X_HAT) - vx.coLinearRapidity() << endl;
    problems = true;
  }


  cout << "{ Tests of boosts from CM frame }\n";

  LorentzVector::setTolerance(.00000000001);
  cout << "    ... tolerance set to " << LorentzVector::getTolerance() << endl;

  v1.set ( Tcomponent(100), -35, -25, 46);
  vx = v1.boostVector();
  v2 = v1.rest4Vector();
  v3 = boostOf(v2,vx);
  if ( v1.isNear(v3) ) {
    cout << "OK: " << v1 << ".rest4Vector = " << v2 <<
	"\n      boosted by " << vx << "-->" << v3 << endl;
  } else {
    cout << "?????" << v1 << ".rest4Vector = " << v2 <<
	"\n      boosted by " << vx << "-->" << v3 << endl;
    problems = true;
  }

  v1.set ( Tcomponent(-100), -23, -25, 46);
  vx = - v1.findBoostToCM();
  v2 = v1.rest4Vector();
  v3 = v2;
  v3.boost(vx);
  if ( v1.isNear(v3) ) {
    cout << "OK: " << v1 << ".rest4Vector = " << v2 <<
	"\n      boosted by " << vx << "-->" << v3 << endl;
  } else {
    cout << "?????" << v1 << ".rest4Vector = " << v2 <<
	"\n      boosted by " << vx << "-->" << v3 << endl;
    problems = true;
  }

  v1.set ( Tcomponent(100), -23, 25, 17);
  vx = - v1.findBoostToCM();
  v2 = v1.rest4Vector();
  v3 = boostOf(v2, 3*vx, vx.beta());
  if ( v1.isNear(v3) ) {
    cout << "OK: " << v1 << ".rest4Vector = " << v2 <<
	"\n      boosted by " << vx << "-->" << v3 << endl;
  } else {
    cout << "?????" << v1 << ".rest4Vector = " << v2 <<
	"\n      boosted by " << vx << "-->" << v3 << endl;
    problems = true;
  }


  cout << "{ Check that Z,Y,X axes give consistent results }\n";

  v1.set ( Tcomponent(100), 23, -25, 37);
  vx.set ( .24, -.31, .4);
  v3 = boostOf(v1,vx);
  cout << "    ... " << v1 << " boosted by " << vx <<
	"\n                    =" << v3 << endl;

  v1.set ( Tcomponent(100), -25, 37, 23);
  vx.set (  -.31, .4, .24);
  {Float8 temp=v3.x(); v3.setX( v3.y() ); v3.setY( v3.z() ); v3.setZ( temp );}
  if ( boostOf(v1,vx).isNear(v3) ) {
    cout << "OK: " << v1 << " boosted by " << vx <<
	"\n                    =" << v3<< endl;
  } else {
    cout << "?????"<< v1 << " boosted by " << vx <<
	"\n                    =" << v3<< endl;
    problems = true;
  }

  v1.set ( Tcomponent(100), 37, 23, -25);
  vx.set (  .4, .24, -.31);
  {Float8 temp=v3.x(); v3.setX( v3.y() ); v3.setY( v3.z() ); v3.setZ( temp );}
  if ( boostOf(v1,vx).isNear(v3) ) {
    cout << "OK: " << v1 << " boosted by " << vx <<
	"\n                    =" << v3<< endl;
  } else {
    cout << "?????"<< v1 << " boosted by " << vx <<
	"\n                    =" << v3<< endl;
    problems = true;
  }


  cout << "{ General boosts that happen to lie on axes }\n";

  LorentzVector::setTolerance(.000000000001);
  cout << "    ... tolerance set to " << LorentzVector::getTolerance() << endl;

  v1.set ( Tcomponent(110), 36, -23, -21);
  vx = .45*X_HAT;
  if ( boostOf(v1,vx).isNear(boostXOf(v1,vx.mag())) ) {
    cout << "OK: " << v1 << " boosted by " << vx <<
	"\n                    =" << boostOf(v1,vx) << endl;
  } else {
    cout << "?????" << v1 << " boosted by " << vx <<
	"\n =" << boostOf(v1,vx) << "!="
	<< boostXOf(v1,vx.mag()) << endl;
    problems = true;
  }


  v1.set ( Tcomponent(120), 34, 23, -21);
  vx = -.35*Y_HAT;
  if ( boostOf(v1,vx).isNear(boostYOf(v1,-vx.mag())) ) {
    cout << "OK: " << v1 << " boosted by " << vx <<
	"\n                    =" << boostOf(v1,vx) << endl;
  } else {
    cout << "?????" << v1 << " boosted by " << vx <<
	"\n =" << boostOf(v1,vx) << "!="
	<< boostYOf(v1,-vx.mag()) << endl;
    problems = true;
  }


  v1.set ( Tcomponent(130), 36, 25, 31);
  vx = .42*Z_HAT;
  if ( boostOf(v1,vx).isNear(boostZOf(v1,vx.mag())) ) {
    cout << "OK: " << v1 << " boosted by " << vx <<
	"\n                    =" << boostOf(v1,vx) << endl;
  } else {
    cout << "?????" << v1 << " boosted by " << vx <<
	"\n =" << boostOf(v1,vx) << "!="
	<< boostZOf(v1,vx.mag()) << endl;
    problems = true;
  }

  cout << "{ Undoing a boost }\n";

  v1.set (10, 20, 30, Tcomponent(40));
  cout << "    ... start from" << v1 << endl;
  vx.set (.3, .4, -.45);
  v2 = boostOf(v1,vx);
  cout << "    ... boosted by " << vx << "-->" << v2 << endl;
  v2.boost(-vx);
  if ( v2.isNear(v1) ) {
    cout << "OK: boost " << -vx << "-->" << v2 << endl;
  } else {
    cout << "????? boost " << -vx << "-->" << v2 << endl;
    problems = true;
  }

} // end of boosts()


void units() {

  cout << "\n --- Unit 4-Vectors\n\n";

  if ( X_HAT4 == LorentzVector ( X_HAT, 0) ) {
    cout << "OK: X_HAT4 = " << X_HAT4 << endl;
  } else {
    cout << "????? X_HAT4 = " << X_HAT4 << endl;
    problems = true;
  }


  if ( Y_HAT4 == LorentzVector ( Y_HAT, 0) ) {
    cout << "OK: Y_HAT4 = " << Y_HAT4 << endl;
  } else {
    cout << "????? Y_HAT4 = " << Y_HAT4 << endl;
    problems = true;
  }


  if ( Z_HAT4 == LorentzVector ( Z_HAT, 0) ) {
    cout << "OK: Z_HAT4 = " << Z_HAT4 << endl;
  } else {
    cout << "????? Z_HAT4 = " << Z_HAT4 << endl;
    problems = true;
  }

  if ( T_HAT4 == LorentzVector ( Tcomponent(1) ) ) {
    cout << "OK: T_HAT4 = " << T_HAT4 << endl;
  } else {
    cout << "????? T_HAT4 = " << T_HAT4 << endl;
    problems = true;
  }

} // end of units()

void input() {

// Input

  cout << "\n --- Input\n\n";

	cout << " { Input of entire vector }\n";

  cin >> v1;
  if (v1.isNear(LorentzVector(99., -99., 95., Tcomponent(500) ))) {
    cout << "OK: input vector is " << v1 << endl;
  } else {
    cout << "????? input vector should be (99, -99, 95; 500) but is " <<
	"\n           " << v1 << endl;
    problems=true;
  }

#ifdef ZAP4
  cin >> v1c;
  // Should fail, and does.
  cout << "New value of const vector v1c is " << v1c << endl;
#endif

	cout << " { Input of cartesian components }\n";

  double xtmp, ytmp, ztmp, ttmp;
  cin >> xtmp >> ytmp >> ztmp >> ttmp ;
  v1.setX(xtmp); v1.setY(ytmp); v1.setZ(ztmp); v1.setT(ttmp);
  if (v1 == LorentzVector(50., 60., 70., Tcomponent(80) )) {
    cout << "OK: input vector is " << v1 << endl;
  } else {
    cout << "????? vector should be ((50, 60, 70), 80) but is "
	"\n                  " << v1 << endl;
    problems=true;
  }

#ifdef ZAP4
  cin >> v1c.x(); // Should fail. Does fail with:
	// svtest.cc:1152: call of overloaded method `operator
	// >>(LorentzVector::WontAlter_X_of_a_const)' is ambiguous
	// candidates are: istream::operator >>(char &)
	// ...
	// operator >>(istream &,LorentzVector::Cin_PutTo_constXPreventer)
#endif

} // end of input()


bool checkAxialRotation ( const LorentzVector &x,
			  const SpaceVector &d,
			  const double& phi,
			  const LorentzVector& result ) {

  // Rotate x by phi around d; return true if result does not match answer
  // within .0001 tolerance; print out appropriate messages.

  LorentzVector xx (x);
  xx.rotate(d, phi);
  if ( xx.isNear(result, .0001) ) {
    cout << "OK: " << x << "Rotated " << phi << " about" << d <<
	"\n = " << xx << "~" << result << endl;
  } else {
    cout << "?????" << x << "Rotated " << phi << " about" << d <<
	"\n = " << xx << " not near " << result << endl;
    return true;
  }

  // Now form a new LorentzVector as xx rotated about the axis by -phi;
  // again check that the result is proper (matches x) very tightly

  LorentzVector yy = rotationOf(xx, d, -phi);
  if ( yy.isNear(x, .000001) ) {
    cout << "OK: " << xx << "Rotated " << -phi << " about" << d <<
	"\n = " << yy << "~" << x << endl;
  } else {
    cout << "?????" << xx << "Rotated " << -phi << " about" << d <<
	"\n = " << yy << " not near " << x << endl;
    return true;
  }

  return false;

  // How did we get those results that we have supplied?
  //
  // We rotate the vector in the following manner:
  // (1) From a set of three coordinate-like vectors with d used as the third:
  //    a - choose V1 = d cross x
  //    b - choose v2 = d cross V1
  //	c - the third is just d
  // (2) Decompose x into u V1 + v V2 + w d
  //	You know a priori that u = 0
  //	Unless d is purely in the Z direction, the two equations dealing with
  // 	the x and y directions suffice:
  //	v = (x.x*d.y - d.x*x.y) / ( V2.x*d.y - d.x*V2.y ) and similar for w.
  // (3) express x in terms of V2-hat instead of V2
  // (4) Rotate in V2-hat, V1hat space by phi
  // (5) rescale back to V1, V2
  // (6) transform back using V1, V2, and d to the x,y,z coordinates.
  // Steps 3--6 combine:
  // x' = - v sin(phi) (|V2| \ |V1|) V1 + v cos(phi) V2 + w d
  //
  // That this result matches the formula used in SpaceVector.cc is a highly
  // non-trivial test.

} /* checkAxialRotation */