// vectest.cc -- general test of all features of PhysVec.h convenience header
//---------------------------------------------------------------------------

#include "ZMutility/ZMenvironment.h"
#include "PhysicsVectors/Rotation.h"
#include "PhysicsVectors/SpaceVector.h"
#include "PhysicsVectors/UnitVector.h"
#include "PhysicsVectors/LorentzVector.h"
#include "PhysicsVectors/LorentzTransformation.h"
#include "PhysicsVectors/PhysVec.h"
ZM_USING_NAMESPACE( zmpv )

USING( std::cout )
USING( std::endl )

int problems( 0 );

  void vv();
  void rv();
  void trans();
  void vr();
  void xyzt();
  void unit();

#define CHECK(CONDITION) 						\
 if ( CONDITION  ) { 							\
   std::cout << "OK: " #CONDITION "\n";					\
 } else {								\
   std::cout << "\n??? " #CONDITION  " >>>>>>>> Line " << __LINE__ << "\n\n";\
   problems++;								\
 }


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 );
}


int main() {

  cout << "\n\n";
  cout << "   *********** \n";
  cout << "   PhysVec\n";
  cout << "   *********** \n";

  vv();
  rv();
  vr();
  xyzt();
  trans();
  unit();

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

} /* End of svtest main() */


void vv() {

  cout << "\n v*v \n\n";

  SpaceVector v1 (2, 3, 4);
  SpaceVector v2 (5, 3, 6);
  CHECK ( v1*v2 == v1.dot(v2) )

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

  CHECK ( w1*w2 == w1.dot(w2) )

}

void rv() {

  cout << "\n r*v \n\n";

  SpaceVector v1 (2, 3, 4);
  Rotation R(.5, .3, .6);
  CHECK ( R*v1 == R(v1) )

  UnitVector u1 (2, -3, 1);
  CHECK ( R*u1 == R(u1) )

  LorentzVector w1 (2, 3, 4, Tcomponent(6));

  if( !(R*w1).isNear(R(w1)) ) {
    cout << "???? R*w1 is NOT near R(w1)" << endl;
    cout << "\nR*w1  = " << R*w1  << endl;
    cout << "\nR(w1) = " << R(w1) << endl;
  } else {
    cout << "OK: R*w1 is near R(w1)" << endl;
  }

  LorentzBoost B ( .2, .3, .4 );
  LorentzTransformation T ( R, B );

  if( !(B*w1).isNear(B(w1)) ) {
    cout << "???? B*w1 is NOT near B(w1)" << endl;
    cout << "\nB*w1  = " << B*w1  << endl;
    cout << "\nB(w1) = " << B(w1) << endl;
  } else {
    cout << "OK: B*w1 is near B(w1)" << endl;
  }

  if( !(T*w1).isNear(T(w1)) ) {
    cout << "???? T*w1 is NOT near T(w1)" << endl;
    cout << "\nT*w1  = " << T*w1  << endl;
    cout << "\nT(w1) = " << T(w1) << endl;
  } else {
    cout << "OK: T*w1 is near T(w1) " << endl;
  }

  RotationX Rx (.4);
  RotationY Ry (.3);
  RotationZ Rz (.2);

  CHECK (Rx*v1 == Rx(v1))
  CHECK (Rx*u1 == Rx(u1))
  CHECK (Rx*w1 == Rx(w1))
  CHECK (Ry*v1 == Ry(v1))
  CHECK (Ry*u1 == Ry(u1))
  CHECK (Ry*w1 == Ry(w1))
  CHECK (Rz*v1 == Rz(v1))
  CHECK (Rz*u1 == Rz(u1))
  CHECK (Rz*w1 == Rz(w1))

  LorentzBoostX Bx (.5);
  LorentzBoostY By (.4);
  LorentzBoostZ Bz (.3);

  CHECK (Bx*w1 == Bx(w1))
  CHECK (By*w1 == By(w1))
  CHECK (Bz*w1 == Bz(w1))

}

void trans() {

  cout << "\n transposeOf()  \n\n";

  Rotation R(.2, .3, .4);
  Rotation RT = inverseOf(R);
  cout << "transposeOf" << R << "is" << RT;
  CHECK ( RT.xx()==R.xx() && RT.xy()==R.yx() && RT.xz()==R.zx() )
  CHECK ( RT.yx()==R.xy() && RT.yy()==R.yy() && RT.yz()==R.zy() )
  CHECK ( RT.zx()==R.xz() && RT.zy()==R.yz() && RT.zz()==R.zz() )

  LorentzBoost B ( .25, .33, .4 );
  LorentzTransformation T ( R, B );
  LorentzTransformation TT = inverseOf(T);
  cout << "contraTransposeOf" << T << "is" << TT << endl;
  CHECK ( TT.xx()==T.xx() && TT.xy()==T.yx() && TT.xz()==T.zx() )
  CHECK ( TT.yx()==T.xy() && TT.yy()==T.yy() && TT.yz()==T.zy() )
  CHECK ( TT.zx()==T.xz() && TT.zy()==T.yz() && TT.zz()==T.zz() )
  CHECK ( TT.xt()== -T.tx() && TT.yt()== -T.ty() )
  CHECK ( TT.zt()== -T.tz() && TT.tt()==  T.tt() )
  CHECK ( T.xt()== -TT.tx() && T.yt()== -TT.ty() )
  CHECK ( T.zt()== -TT.tz() && T.tt()==  TT.tt() )

  TT =  inverseOf(T);
  cout << "transposeOf" << T << "is" << TT << endl;
  CHECK ( TT.xx()==T.xx() && TT.xy()==T.yx() && TT.xz()==T.zx() )
  CHECK ( TT.yx()==T.xy() && TT.yy()==T.yy() && TT.yz()==T.zy() )
  CHECK ( TT.zx()==T.xz() && TT.zy()==T.yz() && TT.zz()==T.zz() )
  CHECK ( TT.xt()==-T.tx() && TT.yt()==-T.ty() )
  CHECK ( TT.zt()==-T.tz() && TT.tt()== T.tt() )
  CHECK ( T.xt()==-TT.tx() && T.yt()==-TT.ty() )
  CHECK ( T.zt()==-TT.tz() && T.tt()== TT.tt() )

}

// Earlier versions had defined v * R as R.transpose  v.
//
// This is so deceptive -- as is the possible alternative interpretation --
// that we simply removed the operation.

void vr() {

  cout << "\n v*r \n\n";

  SpaceVector v1 (2, 3, 4);
  SpaceVector v2 (1, 5, -4);
  Rotation R(.5, .3, .6);
  CHECK ( (v1*R).isNear((inverseOf(R))(v1)) )
  CHECK ( close((v1*R)*v2, v1.dot(R(v2))) )

  UnitVector u1 (2, -3, 1);
  UnitVector u2 (-2, -2, 4);
  CHECK ( (u1*R).isNear((inverseOf(R))(u1)) )
  CHECK ( close((u1*R)*u2, u1.dot(R(u2))) )

  LorentzVector w1 (2, 3, 4, Tcomponent(6));
  LorentzVector w2 (-2, 3, -3, Tcomponent(8));
  CHECK ( (w1*R).isNear((inverseOf(R))(w1)) )
  CHECK ( close((w1*R)*w2, w1.dot(R(w2))) )

  LorentzBoost B ( .2, .3, .4 );
  LorentzTransformation T ( R, B );
  CHECK ( (w1*B).isNear((inverseOf(B))(w1)) )
  CHECK ( close((w1*B)*w2, w1.dot(B(w2))) )
  CHECK ( (w1*T).isNear((inverseOf(T))(w1)) )
  CHECK ( close((w1*T)*w2, w1.dot(T(w2))) )

  RotationX Rx (.4);
  RotationY Ry (.3);
  RotationZ Rz (.2);

  CHECK ( (v1*Rx).isNear((inverseOf(Rx))(v1)) )
  CHECK ( close((v1*Rx)*v2, v1.dot(Rx(v2))) )
  CHECK ( (v1*Ry).isNear((inverseOf(Ry))(v1)) )
  CHECK ( close((v1*Ry)*v2, v1.dot(Ry(v2))) )
  CHECK ( (v1*Rz).isNear((inverseOf(Rz))(v1)) )
  CHECK ( close((v1*Rz)*v2, v1.dot(Rz(v2))) )

  CHECK ( (u1*Rx).isNear((inverseOf(Rx))(u1)) )
  CHECK ( close((u1*Rx)*u2, u1.dot(Rx(u2))) )
  CHECK ( (u1*Ry).isNear((inverseOf(Ry))(u1)) )
  CHECK ( close((u1*Ry)*u2, u1.dot(Ry(u2))) )
  CHECK ( (u1*Rz).isNear((inverseOf(Rz))(u1)) )
  CHECK ( close((u1*Rz)*u2, u1.dot(Rz(u2))) )

  CHECK ( (w1*Rx).isNear((inverseOf(Rx))(w1)) )
  CHECK ( close((w1*Rx)*w2, w1.dot(Rx(w2))) )
  CHECK ( (w1*Ry).isNear((inverseOf(Ry))(w1)) )
  CHECK ( close((w1*Ry)*w2, w1.dot(Ry(w2))) )
  CHECK ( (w1*Rz).isNear((inverseOf(Rz))(w1)) )
  CHECK ( close((w1*Rz)*w2, w1.dot(Rz(w2))) )

  LorentzBoostX Bx (.5);
  LorentzBoostY By (.4);
  LorentzBoostZ Bz (.3);

  #ifdef EXCISED // removed after CLHEP merge. Dont even know why w1*Rx works...
  CHECK ( (w1*Bx).isNear((inverseOf(Bx))(w1)) )
  CHECK ( close((w1*Bx)*w2, w1.dot(Bx(w2))) )
  CHECK ( (w1*By).isNear((inverseOf(By))(w1)) )
  CHECK ( close((w1*By)*w2, w1.dot(By(w2))) )
  CHECK ( (w1*Bz).isNear((inverseOf(Bz))(w1)) )
  CHECK ( close((w1*Bz)*w2, w1.dot(Bz(w2))) )
  #endif

}

void xyzt() {

  cout << "\n v*r \n\n";

}

void unit() {

}