#include <iostream>
using std::cout;
using std::cerr;
using std::cin;
using std::endl;

#include <valarray>
using std::valarray;

#include <fstream>
using std::fstream;

#include <string>
using std::string;

#include <list>
using std::list;

#include <vector>
using std::vector;

#include "TObject.h"
#include "TVector3.h"
#include "TRandom.h"
#include "TMatrixD.h"
#include "TVectorD.h"
#include "TDecompSVD.h"
#include "Math/Vector3D.h"
using ROOT::Math::XYZVector;

#include "Math/Point3D.h"
using ROOT::Math::XYZPoint;

#include "Math/Transform3D.h"
using ROOT::Math::Transform3D;

#include "Math/RotationX.h"
using ROOT::Math::RotationX;
#include "Math/RotationY.h"
using ROOT::Math::RotationY;
#include "Math/RotationZ.h"
using ROOT::Math::RotationZ;
#include "Math/Rotation3D.h"
using ROOT::Math::Rotation3D;

#include "PndDrcPhoton.h"
#include "PndDrcOptReflAbs.h"
#include "PndDrcSurfAbs.h"
#include "PndDrcSurfPolyFlat.h"
#include "PndDrcSurfQuadFlatDiff.h"
#include "PndDrcOptMatAbs.h"
#include "PndDrcOptDev.h"
#include "PndDrcUtil.h"

//----------------------------------------------------------------------
PndDrcSurfQuadFlatDiff::PndDrcSurfQuadFlatDiff()
{
  m_s1 = 0;
  m_s2 = 0;
}
//----------------------------------------------------------------------
PndDrcSurfQuadFlatDiff::~PndDrcSurfQuadFlatDiff()
{
  if (m_s1) delete m_s1;
  if (m_s2) delete m_s2;
}

//----------------------------------------------------------------------
PndDrcSurfQuadFlatDiff* PndDrcSurfQuadFlatDiff::clone() const
{
  return new PndDrcSurfQuadFlatDiff(*this); 
}
//----------------------------------------------------------------------
void PndDrcSurfQuadFlatDiff::copy(const PndDrcSurfQuadFlatDiff& s)
{
  m_s1p1        = s.m_s1p1;
  m_s1p2        = s.m_s1p2;
  m_s2p1        = s.m_s2p1;
  m_s2p2        = s.m_s2p2;
  m_s1          = (s.m_s1)->clone();
  m_s2          = (s.m_s2)->clone();
  m_normal      = s.m_normal;
  m_surfAux     = s.m_surfAux;
}
//----------------------------------------------------------------------
PndDrcSurfQuadFlatDiff::PndDrcSurfQuadFlatDiff(const PndDrcSurfQuadFlatDiff& s) : PndDrcSurfAbs(s)
{
  if (s.m_verbosity>=1) cout<<"  PndDrcSurfQuadFlatDiff::PndDrcSurfQuadFlatDiff"
			    <<"(const PndDrcSurfQuadFlatDiff&) name,copy: "  
			    <<s.m_name<<" "<<s.m_copyNumber<<endl;
  copy(s);
}
//----------------------------------------------------------------------
PndDrcSurfQuadFlatDiff&  PndDrcSurfQuadFlatDiff::operator=(const PndDrcSurfQuadFlatDiff& s)
{
  if (s.m_verbosity>=1) cout<<"  PndDrcSurfQuadFlatDiff::operator="
			    <<"(const PndDrcSurfQuadFlatDiff&) name,copy: "  
			    <<s.m_name<<" "<<s.m_copyNumber<<endl;
  if (&s != this)
    {
      static_cast<PndDrcSurfAbs&>((*this)) = s; // assignment of base class part.
      copy(s);
    }

  return *this;
}
//----------------------------------------------------------------------
void PndDrcSurfQuadFlatDiff::addSurface(const PndDrcSurfAbs& surf, XYZPoint p1, XYZPoint p2)
{
  if (m_s1 && m_s2)
    {
      cerr<<" *** PndDrcSurfQuadFlatDiff::addSurface: you can only add two surfaces.\n";
      exit(EXIT_FAILURE);
    }
  else if (!m_s1)  
    {
      m_s1 = surf.clone();
      m_s1p1 = p1;
      m_s1p2 = p2;
    }
  else if (!m_s2)
    {
      m_s2 = surf.clone();
      m_s2p1 = p1;
      m_s2p2 = p2;
    }
  else
    {
      cerr<<" *** PndDrcSurfQuadFlatDiff::addSurface: this line should never hit."<<endl;
      exit(EXIT_FAILURE);
    }
  
  if (m_s1 && m_s2)
    {
      XYZVector cross1 = (m_s1p2-m_s1p1).Cross(m_s2p1-m_s1p1);
      XYZVector cross2 = (m_s1p2-m_s1p1).Cross(m_s2p2-m_s1p1);
      m_normal = (cross1.Mag2() > cross2.Mag2()) ? cross1.Unit() : cross2.Unit();
      m_surfAux.addPoint(m_s1p1);
      m_surfAux.addPoint(m_s1p2);
      m_surfAux.addPoint(m_s2p1);
      m_surfAux.addPoint(m_s2p2);
    }
}
//----------------------------------------------------------------------
XYZVector PndDrcSurfQuadFlatDiff::normal(const XYZPoint& point) const
{
  // normal vector is a constant.
  return m_normal;
}
//----------------------------------------------------------------------
bool PndDrcSurfQuadFlatDiff::surfaceHit(PndDrcPhoton& ph, 
				     XYZPoint&  pos_new, 
				     double&    path_length) const
{
  static const double eps = 1.0e-9;

  if (verbosity()>=3) cout<<"      PndDrcSurfQuadFlatDiff::surfaceHit: name="
			  <<name()<<endl;


  if (!m_s1 || !m_s2)
    {
      cerr<<" *** PndDrcSurfQuadFlatDiff::surfaceHit: at least one of"
	  <<" both surfaces undefined."<<endl;
      exit(EXIT_FAILURE);
    }

 

  // find first point

  XYZPoint  point = ph.position();
  XYZVector dir   = ph.direction().Unit();


  // calculate hits with area
  // (point+lambda*dir - m_p[0]) * m_normal = 0
  //
  // point*m_normal + lambda*dir*m_normal - m_p[0]*m_normal = 0
  //
  // lamda = (m_p[0]-point)*m_normal  / dir*m_normal
  //
  // dir*m_normal > 0 , see condition above

  double h = dir.Dot(m_normal);
  if (h==0) 
    {
      if (verbosity()>=4) cout<<"      PndDrcSurfQuadFlatDiff::surfaceHit:"
			      <<" dir perp norm, no hit. \n";
     return false;
    }
  double   lambda = (m_s1p1-point).Dot(m_normal) / dir.Dot(m_normal);
  pos_new = point + lambda*dir;

  if (verbosity()>=4) cout<<"      PndDrcSurfQuadFlatDiff::surfaceHit:"<<name()<<endl
			  <<"      lambda pos_new ="<<lambda<<" "<<pos_new<<endl;


  if (lambda<-eps) // flying in wrong direction.
    {
      if (verbosity()>=4) cout<<"      PndDrcSurfQuadFlatDiff::surfaceHit: wrong dir "
			      <<" lambda par ="<<lambda<<" no hit."<<endl;
      return false; 
    }

  path_length = fabs(lambda);

  // check surfaces between rectangular partial surface and eventual bended
  // surface borders.
  bool hit_aux_surf = m_surfAux.withinSurface(pos_new);

  // if both are flat, there was not hit.
  if (m_s1->isFlat() && m_s2->isFlat()) 
    {
      if (hit_aux_surf)
	{
	  if (m_pixel) ph.setFate(Drc::PhotMeasured);
	    
	  if (verbosity()>=4) cout<<"      PndDrcSurfQuadFlatDiff::surfaceHit:"
				  <<" within rectangular partial surface, hit. "
				  <<endl;
	  return true;
	}
      else
	{
	  if (verbosity()>=4) cout<<"      PndDrcSurfQuadFlatDiff::surfaceHit:"
				  <<" no hit within reactangular partial surface"
				  <<" of flat borders. "
				  <<endl;
	  return false;
	}
    }


  bool in_between1 = false;
  if (!(m_s1->isFlat()))
    {
      // check if point between curve and line
      // create a fake photon.
      if (verbosity()>=4) cout<<"      PndDrcSurfQuadFlatDiff::surfaceHit: (1) "
			      <<" check point between line and curve of "
			      <<m_s1->name()<<endl;
      PndDrcPhoton ph1;
      ph1.setPosition(pos_new);

      XYZPoint ref(0.5*(XYZVector(m_s1p1) + XYZVector(m_s1p2)));
      //                 no addition between XYZPoints possible
      ref.SetX(1.01*ref.X()); // bend it a little bit outside
      ref.SetY(1.01*ref.Y()); // so that dir points a little bit inside.
      XYZVector dir = (XYZPoint(pos_new)-ref).Unit();

      ph1.setDirection(dir);

     if (verbosity()>=4) 
	{
	  cout<<"        with aux. photon "<<endl;
	  cout<<"        ph position "<<ph1.position()<<endl;
	  cout<<"        & direction "<<ph1.direction()<<endl;
	}

      XYZPoint v_dummy;
      double dummy;
      if (m_s1->surfaceHit(ph1,v_dummy,dummy)) 
	{
	  if (verbosity()>=4) cout<<"      PndDrcSurfQuadFlatDiff::surfaceHit: (1) "
				  <<" between line and curve of "
				  <<m_s1->name()<<", hit."<<endl;
	  in_between1 =  true;
	}
      else
	{
	  if (verbosity()>=4) cout<<"      PndDrcSurfQuadFlatDiff::surfaceHit: "
				  <<" no hit between line and curve of "
				  <<m_s1->name()<<", check other side."<<endl;
	}
    }


  //                   convex           concave  cases
  bool in_between2 = (hit_aux_surf && m_s2->isFlat()); // for later logic
  if (! m_s2->isFlat())
    {
      // check if point between curve and line
      // create a fake photon.
      if (verbosity()>=4) cout<<"      PndDrcSurfQuadFlatDiff::surfaceHit: (2) "
			      <<" check point between line and curve of "
			      <<m_s2->name()<<endl;
      PndDrcPhoton ph1;
      ph1.setPosition(pos_new);

      XYZPoint ref(0.5*(XYZVector(m_s2p1) + XYZVector(m_s2p2)));
      //                 no addition between XYZPoints possible
      ref.SetX(1.01*ref.X()); // bend it a little bit outside
      ref.SetY(1.01*ref.Y()); // so that dir points a little bit inside.
      XYZVector dir = (XYZPoint(pos_new)-ref).Unit();

      ph1.setDirection(dir);
 
     if (verbosity()>=4) 
	{
	  cout<<"        with aux. photon "<<endl;
	  cout<<"        ph position "<<ph1.position()<<endl;
	  cout<<"        & direction "<<ph1.direction()<<endl;
	}

      XYZPoint v_dummy;
      double dummy;
      // m_s2->setVerbosity(4);//###
      if (m_s2->surfaceHit(ph1,v_dummy,dummy)) 
	{
	  if (verbosity()>=4) cout<<"      PndDrcSurfQuadFlatDiff::surfaceHit: (2) "
				  <<" between line and curve of "
				  <<m_s2->name()<<", hit."<<endl;
	  in_between2 = true;
	}
      else
	{
	  if (verbosity()>=4) cout<<"      PndDrcSurfQuadFlatDiff::surfaceHit: "
				  <<" no hit between line and curve of "
				  <<m_s2->name()<<", no hit possibility left."<<endl;
	}
    }


  if (verbosity()>=4) 
    {
      cout<<"      PndDrcSurfQuadFlatDiff::surfaceHit: "<<endl
	  <<"         hit_aux_surf = "<<hit_aux_surf<<endl;
      if (m_s1->isFlat()) cout<<"         in_between1  = "<<"flat"<<endl;
      if (!m_s1->isFlat()) cout<<"         in_between1  = "<<in_between1<<endl;
      if (m_s2->isFlat()) cout<<"         in_between2  = "<<"flat"<<endl;
      if (!m_s2->isFlat()) cout<<"         in_between2  = "<<in_between2<<endl;

    }


  // in/exclusion logic
  if (hit_aux_surf)
    {
      // exclude photons from concave parts
      if ( !(m_s1->isFlat()) && in_between1) return false;
      if ( !(m_s2->isFlat()) && in_between2) return false;
      return true;
    }
  else
    {
      if (in_between1) return true; // surface convex, photon inside
      if (in_between2) return true; // surface convex, photon inside
      return false;
    }
}
//----------------------------------------------------------------------
void PndDrcSurfQuadFlatDiff::print(fstream& stream) const
{
  stream<<" TPolyLine3D *l = new TPolyLine3D(2);"<<endl;

  stream<<" l->SetPoint("<<0<<","
	<<m_s1p1.X()<<","
	<<m_s1p1.Y()<<","
	<<m_s1p1.Z()<<");"<<endl;
  stream<<" l->SetPoint("<<1<<","
	<<m_s2p1.X()<<","
	<<m_s2p1.Y()<<","
	<<m_s2p1.Z()<<");"<<endl;
  stream<<" l->SetLineColor("<<m_printColor<<");"<<endl;
  stream<<" l->Draw();"<<endl;

  stream<<" TPolyLine3D *l = new TPolyLine3D(2);"<<endl;

  stream<<" l->SetPoint("<<0<<","
	<<m_s1p2.X()<<","
	<<m_s1p2.Y()<<","
	<<m_s1p2.Z()<<");"<<endl;
  stream<<" l->SetPoint("<<1<<","
	<<m_s2p2.X()<<","
	<<m_s2p2.Y()<<","
	<<m_s2p2.Z()<<");"<<endl;
  stream<<" l->SetLineColor("<<m_printColor<<");"<<endl;
  stream<<" l->Draw();"<<endl;
}
//----------------------------------------------------------------------
void PndDrcSurfQuadFlatDiff::addTransform(const Transform3D& trans)
{
  if (!m_s1 || !m_s2) 
    {
      cerr<<" PndDrcSurfQuadFlatDiff::shift: surface undefined"<<endl;
      exit(EXIT_FAILURE);
    }

  if (verbosity()>=3) cout<<" PndDrcSurfQuadFlatDiff::addTransform() name="<<m_name<<endl;

  m_normal = trans * m_normal;
  m_s1p1   = trans * m_s1p1;
  m_s1p2   = trans * m_s1p2;
  m_s2p1   = trans * m_s2p1;
  m_s2p2   = trans * m_s2p2;
  m_s1->addTransform(trans);
  m_s2->addTransform(trans);
  m_surfAux.addTransform(trans);
}