// SurfDCA_t.cpp

// test SurfDCA.

#include "SurfDCA.h"
#include <assert.h>
#include <bstring.h>
#include <sstream>
#include <iostream>
#include "objstream/StdObjStream.hpp"
#include "trfutil/trfstream.h"
#include "trfutil/TRFMath.h"
#include "spacegeom/SpacePath.h"
#include "trfbase/CrossStat.h"
#include "trfbase/VTrack.h"
#include "trfbase/SurfTest.h"

using std::cout;
using std::cerr;
using std::endl;
using std::istringstream;
using std::ostringstream;

using namespace trf;

//**********************************************************************
namespace{
// comparison of doubles
int myequal(double x1, double x2) {
  double small = 1.e-12;
  if ( fabs(x1-x2) < small ) return 1;
  cout << "myequal: difference too large:" << endl;
  cout << "value 1: " << x1 << endl;
  cout << "value 2: " << x2 << endl;
  cout << "   diff: " << x1-x2 << endl;
  cout << "maxdiff: " << small << endl;
  return 0;
}
}
//**********************************************************************

int main( ) {

  string ok_prefix = "test SurfDCA (I): ";
  string error_prefix = "test SurfDCA (E): ";

  cout << ok_prefix
       << "------ Testing component SurfDCA. ------" << endl;

  // Make sure assert is enabled.
  bool assert_flag = false;
  assert ( ( assert_flag = true, assert_flag ) );
  if ( ! assert_flag ) {
    cerr << "Assert is disabled" << endl;
    return 1;
  }

  //********************************************************************

  ObjTable::register_type<SurfDCA>();

  cout << trf_format;
  cout << ok_prefix << "Check indices." << endl;
  assert( SurfDCA::IRSIGNED == 0 );
  assert( SurfDCA::IZ       == 1 );
  assert( SurfDCA::IPHID    == 2 );
  assert( SurfDCA::ITLM     == 3 );
  assert( SurfDCA::IQPT     == 4 );

  //********************************************************************

  cout << ok_prefix << "Test constructor." << endl;
  SurfDCA sdca1;
  cout << sdca1 << endl;
  SurfDCA sdca_x1(0.1,0.2);
  SurfDCA sdca_x2(0.1,0.3);
  cout << sdca_x1 << endl;

  //********************************************************************

  cout << ok_prefix << "Test equality and ordering." << endl;
  assert( sdca1.pure_equal(sdca1) );
  assert( ! sdca1.pure_less_than(sdca1) );
  assert( sdca_x1.pure_equal(sdca_x1) );
  assert( ! sdca_x1.pure_less_than(sdca_x1) );
  assert( !sdca_x1.pure_equal(sdca1) );
  assert( sdca1.pure_less_than(sdca_x1) );
  assert( sdca_x1.pure_less_than(sdca_x2) );

  //********************************************************************

  cout << ok_prefix << "Test virtual constructor." << endl;
  SurfDCA* psdca4 = sdca1.new_pure_surface();
  if ( sdca1 != *psdca4 ) {
    cerr << error_prefix << "Virtual construction failed." << endl;
    return 5;
  }
  SurfDCA* psdca_x4 = sdca_x1.new_pure_surface();
  if ( sdca_x1 != *psdca_x4 ) {
    cerr << error_prefix << "Virtual construction failed." << endl;
    return 5;
  }

  //********************************************************************

  cout << ok_prefix << "Test type." << endl;
  cout << SurfDCA::get_static_type() << endl;
  // method get_type is implemented in class Surface
  // it returns get_pure_type, which is implemented in SurfDCA 
  cout << sdca1.get_type() << endl;
  assert( sdca1.get_type() != 0 );
  assert( sdca1.get_type() == SurfDCA::get_static_type() );

  //********************************************************************

  cout << ok_prefix << "Test parameter access." << endl;

  assert( sdca_x1.get_parameter(SurfDCA::IX) == 0.1);
  assert( sdca_x1.get_parameter(SurfDCA::IY) == 0.2);
  assert( sdca1.get_parameter(SurfDCA::IX) == 0.0);
  assert( sdca1.get_parameter(SurfDCA::IY) == 0.0);

  //********************************************************************

  cout << ok_prefix << "construct a Track Vector." << endl;

  TrackVector tvec;
  tvec(0) = 1.0;                                // r_signed
  tvec(1) = 2.0;                                // z
  tvec(2) = 3.0;                                // phi_direction
  tvec(3) = 4.0;                                // tlm
  tvec(4) = 5.0;                                // qpt
  cout << " *** TrackVector tvec = " << tvec << endl;
  cout << " *** r_signed      = " << tvec(0) << endl;
  cout << " *** z             = " << tvec(1) << endl;
  cout << " *** phi_direction = " << tvec(2) << endl;
  cout << " *** tlm           = " << tvec(3) << endl;
  cout << " *** qpt           = " << tvec(4) << endl;

  double r    = abs(tvec(0));
  double z    = tvec(1);
  double sign = 0.0;
  double phi  = 0.0;
  if ( tvec(0) != 0.0 ) {
    sign = tvec(0)/abs(tvec(0));
    phi  = tvec(2)-(sign*PI2);
    phi  = fmod2( phi, TWOPI );
  }
  double tlam = tvec(3);

  cout << " *** r   = " << r   << endl;
  cout << " *** phi = " << phi << endl;
  cout << " *** z   = " << z   << endl;

  assert( fabs(phi-( 1.4292)) < 1.e-4 ) ;      // for r_signed = +1
  //  assert( fabs(phi-(-1.7124)) < 1.e-4 ) ;      // for r_signed = -1


  //********************************************************************

  cout << ok_prefix << "Test crossing status." << endl;

  double clam = 1.0/sqrt(1.0+tlam*tlam);
  double slam = tlam/sqrt(1.0+tlam*tlam);

  double dz_ds = slam;
  
  double alpha_in     = -2.0;
  double calf_in      = cos(alpha_in);
  double salf_in      = sin(alpha_in);
  double dr_ds_in     = clam*calf_in;
  double r_dphi_ds_in = clam*salf_in;
  CylindricalPath cpth_in(r, phi, z, dr_ds_in, r_dphi_ds_in, dz_ds);
  SurfTest stest_in(cpth_in);

  double alpha_out     =  1.0;
  double calf_out      = cos(alpha_out);
  double salf_out      = sin(alpha_out);
  double dr_ds_out     = clam*calf_out;
  double r_dphi_ds_out = clam*salf_out;
  CylindricalPath cpth_out(r, phi, z, dr_ds_out, r_dphi_ds_out, dz_ds);
  SurfTest stest_out(cpth_out);

  VTrack ton( SurfacePtr(sdca1.new_pure_surface()), tvec );
  cout << " sdca1       = " << sdca1 << endl;
  cout << " VTrack ton  = " << ton << endl;
  VTrack tin( SurfacePtr(stest_in.new_pure_surface()), tvec );
  cout << " stest_in    = " << stest_in << endl;
  cout << " VTrack tin  = " << tin << endl;
  VTrack tout( SurfacePtr(stest_out.new_pure_surface()), tvec );
  cout << " stest_out   = " << stest_out << endl;
  cout << " VTrack tout = " << tout << endl;

  SpacePath sp_on  = ton.space_vector();
  SpacePath sp_in  = tin.space_vector();
  SpacePath sp_out = tout.space_vector();

  cout << " *** sp_on *** " << endl;
  cout << sp_on  << endl;
  cout << " *** sp_in *** " << endl;
  cout << sp_in  << endl;
  cout << " *** sp_out *** " << endl;
  cout << sp_out << endl;

  CrossStat xs1 = sdca1.pure_status(ton);
  cout << "xs1 = " << xs1 << endl;
  assert( xs1.at() && xs1.on() && !xs1.inside() && !xs1.outside()
          && !xs1.in_bounds() && ! xs1.out_of_bounds() );

  cout << stest_in << endl;
  CrossStat xs2 = sdca1.pure_status(tin);
  cout << "xs2 = " << xs2 << endl;
  assert( !xs2.at() && !xs2.on() && xs2.inside() && !xs2.outside()
          && !xs2.in_bounds() && ! xs2.out_of_bounds() );

  cout << stest_out << endl;
  CrossStat xs3 = sdca1.pure_status(tout);
  cout << "xs3 = " << xs3 << endl;
  assert( !xs3.at() && !xs3.on() && !xs3.inside() && xs3.outside()
          && !xs3.in_bounds() && ! xs3.out_of_bounds() );

  //********************************************************************

  cout << ok_prefix << "Test vector difference." << endl;
  TrackVector tvec2;
  tvec2(0) = 1.1;                                // r_signed
  tvec2(1) = 2.2;                                // z
  tvec2(2) = 3.3;                                // phi_direction
  tvec2(3) = 4.4;                                // tlm
  tvec2(4) = 5.5;                                // qpt
  TrackVector diff = sdca1.vec_diff(tvec2,tvec);
  TrackVector ediff;
  ediff(0) = 0.1;
  ediff(1) = 0.2;
  ediff(2) = 0.3;
  ediff(3) = 0.4;
  ediff(4) = 0.5;
  TrackVector zero = diff - ediff;
  if ( zero.amax() > 1.e-10 )  {
    cerr << error_prefix << "Incorrect difference." << endl;
    return 9;
  }
  //  assert( zero.amax() <= 1.e-10 );

  //********************************************************************

  cout << ok_prefix << "Test space point." << endl;
  SpacePoint spt = ton.space_point();
  assert( spt.rxy() == 1.0 );
  assert( fabs(spt.phi()-phi) < 1.e-4 );
  assert( spt.z()   == 2.0 );

  TrackVector vec_1;
  vec_1(SurfDCA::IZ)=0.3;
  VTrack trv_1(SurfacePtr(new SurfDCA(0.1,0.2)));
  trv_1.set_vector(vec_1);
  trv_1.set_forward();
  assert( myequal(trv_1.space_point().x(),0.1) );
  assert( myequal(trv_1.space_point().y(),0.2) );
  assert( myequal(trv_1.space_point().z(),0.3) );
  vec_1(SurfDCA::IRSIGNED)=0.3;
  vec_1(SurfDCA::IPHID)=PI2;
  trv_1.set_vector(vec_1);
  trv_1.set_forward();
  assert( myequal(trv_1.space_point().x(),0.1+0.3) );
  assert( myequal(trv_1.space_point().y(),0.2) );
  assert( myequal(trv_1.space_point().z(),0.3) );

  vec_1(SurfDCA::IRSIGNED)=0.3;
  vec_1(SurfDCA::IPHID)=-PI;
  trv_1.set_vector(vec_1);
  trv_1.set_forward();
  assert( myequal(trv_1.space_point().x(),0.1) );
  assert( myequal(trv_1.space_point().y(),0.2+0.3) );
  assert( myequal(trv_1.space_point().z(),0.3) );
  

  //********************************************************************

  cout << ok_prefix << "Test space vector." << endl;
  SpacePath svec = ton.space_vector();
  double lam = atan(4.0);
  assert( svec.rxy() == 1.0 );
  assert( fabs(spt.phi()-phi) < 1.e-4 );
  assert( svec.z()   == 2.0 );
  assert( myequal( svec.drxy(),     cos(lam)*cos(PI2) ) );
  assert( myequal( svec.rxy_dphi(), cos(lam)*sin(sign*PI2) ) );
  assert( myequal( svec.dz(),       sin(lam)          ) );

  vec_1(SurfDCA::IRSIGNED)=0.3;
  vec_1(SurfDCA::IPHID)=-PI;
  trv_1.set_vector(vec_1);
  trv_1.set_forward();
  assert( myequal(trv_1.space_vector().x(),0.1) );
  assert( myequal(trv_1.space_vector().y(),0.2+0.3) );
  assert( myequal(trv_1.space_vector().z(),0.3) );
  assert( myequal(trv_1.space_vector().dx(),-1.0) );
  assert( myequal(trv_1.space_vector().dz(),0.0) );
  assert( myequal(trv_1.space_vector().dy(),0.0) );
  vec_1(SurfDCA::IRSIGNED)=0.3;
  vec_1(SurfDCA::IPHID)=PI2;
  trv_1.set_vector(vec_1);
  trv_1.set_forward();
  assert( myequal(trv_1.space_vector().x(),0.1+0.3) );
  assert( myequal(trv_1.space_vector().y(),0.2) );
  assert( myequal(trv_1.space_vector().z(),0.3) );
  assert( myequal(trv_1.space_vector().dx(),0.0) );
  assert( myequal(trv_1.space_vector().dz(),0.0) );
  assert( myequal(trv_1.space_vector().dy(),1.0) );


  //********************************************************************

  cout << ok_prefix << "Write object data." << endl;
  {
    ObjPtr pobj( new SurfDCA );
    ostringstream mystream;
    StdObjStream objstream(mystream);
    objstream.write_object("dca1",pobj);
    cout << mystream.str() << endl;
    assert( ObjTable::has_object_name("dca1") );
    string::size_type pos = 0;
    assert( (pos=mystream.str().find( "dca1 ",pos)) != string::npos );
    assert( (pos=mystream.str().find( "SurfDCA ",pos)) != string::npos );
  }

  //********************************************************************

  cout << ok_prefix << "Read object data." << endl;
  {
    string istring = "[ dca2 SurfDCA x=0.1 y=0.2 ]";
    istringstream isstrm(istring);
    {
      StdObjStream obstr(isstrm);
      string name = obstr.read_object();
      assert( name == "dca2" );
      const SurfDCA* pobj;
      ObjTable::get_object(name,pobj);
      assert( pobj != 0 );
      assert( pobj->get_type() == SurfDCA::get_static_type() );
      assert( pobj->get_x() == 0.1);
      assert( pobj->get_y() == 0.2);
    }
  }

  //********************************************************************

  cout << ok_prefix
       << "------------- All tests passed. -------------" << endl;
  return 0;

  //********************************************************************

}