#include "ThinZPlaneMS.h"
#include "trfbase/ETrack.h"
#include "SurfZPlane.h"
#include "trfbase/SurfacePtr.h"
#include "trfutil/trfstream.h"
#include <cassert>
#include <string>
#include <sstream>
#include "objstream/StdObjStream.hpp"

#ifndef DEFECT_NO_STDLIB_NAMESPACES
using std::cout;
using std::cerr;
using std::endl;
using std::string;
using std::ostringstream;
using std::istringstream;
#endif

using namespace trf;

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

int main( ) {

  string component = "ThinZPlaneMS";
  string ok_prefix = component + " (I): ";
  string error_prefix = component + " test (E): ";

  cout << ok_prefix
       << "---------- Testing component " + component
          + ". ----------" << 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;
  }

  //********************************************************************
  cout << trf_format;

  ThinZPlaneMS scatterIt(1.0);

  TrackError initialError;
  initialError(2,2) = 0.0;
  initialError(3,3) = 0.0;
  initialError(4,4) = 0.0;
  initialError(4,3) = 0.0;
  initialError(3,4) = 0.0;

  TrackVector trv;
  trv(0) = 1.0;
  trv(1) = 1.0;
  trv(2) = 1.0;
  trv(3) = 1.0;
  trv(4) = 1.0;

  SurfacePtr srf(new SurfZPlane( 20.0 ));

  ETrack tmpTrack( srf, trv, initialError );
  const TrackError& finalError = tmpTrack.get_error();

  tmpTrack.set_error( initialError );
  scatterIt.interact( tmpTrack );

  double f = trv(2);
  double g = trv(3);

  double pi = acos(-1.0);
  double theta = atan(sqrt(f*f + g*g));
  double phi = atan(sqrt((g*g)/(f*f)));

  double xhat = sin(theta)*cos(phi);
  double yhat = sin(theta)*sin(phi);
  double zhat = sqrt(1 - sin(theta)*sin(theta));

  double radLength = scatterIt.get_rad_length();

  double trueLength = radLength/cos(theta);

  double thetaSqr = (0.0136)*1.0*sqrt(trueLength)*
                         (1 + 0.038*log(trueLength));
  thetaSqr *= thetaSqr;

  double norm = sqrt(xhat*xhat + yhat*yhat);

  double dxzxz = (yhat/(norm*zhat))*(yhat/(norm*zhat));
  dxzxz += pow((xhat/norm)*(1 + (norm*norm)/(zhat*zhat)),2.0);
  double dyzyz = (xhat/(norm*zhat))*(xhat/(norm*zhat));
  dyzyz += pow((yhat/norm)*(1 + (norm*norm)/(zhat*zhat)),2.0);
  double dxzyz = -xhat*yhat/(zhat*zhat);
  dxzyz += xhat*yhat/(norm*norm)*pow((1 + (norm*norm)/(zhat*zhat)),2.0);
 
  dxzxz *= thetaSqr;
  dyzyz *= thetaSqr;
  dxzyz *= thetaSqr;
 
  assert(scatterIt.get_rad_length() == 1.0);
  assert(finalError(0,0)==0.0);
  assert(finalError(1,1)==0.0);
  assert((finalError(2,2)<(dxzxz+.00001))&&((finalError(2,2)>(dxzxz-.00001))));
  assert((finalError(3,3)<(dyzyz+.00001))&&((finalError(3,3)>(dyzyz-.00001))));
  assert(finalError(4,4)==0.0);
  assert((finalError(2,3)<(dxzyz+.00001))&&((finalError(2,3)>(dxzyz-.00001))));
  assert(finalError(2,3)==finalError(3,2));


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

  ObjTable::register_type<ThinZPlaneMS>();

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

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

  cout << ok_prefix << "Read object data." << endl;
  {
    string istring = "[ obj2 ThinZPlaneMS rad_length=1.11; ]";
    istringstream isstrm(istring);
    {
      StdObjStream obstr(isstrm);
      string name = obstr.read_object();
      assert( name == "obj2" );
      const ThinZPlaneMS* pobj;
      ObjTable::get_object(name,pobj);
      assert( pobj != 0 );
      assert( pobj->get_type() == ThinZPlaneMS::get_static_type() );
      assert( pobj->get_rad_length() == 1.11 );
    }
  }

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

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

  //********************************************************************
 
}