#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 <limits>
using std::numeric_limits;

#include "TObject.h"
#include "TVector3.h"
#include "TRandom.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 "DrcPhoton.h"
#include "DrcOptReflAbs.h"
#include "DrcSurfAbs.h"
#include "DrcOptMatAbs.h"
#include "DrcOptDev.h"
#include "DrcOptVol.h"
#include "DrcUtil.h"

//----------------------------------------------------------------------
DrcOptVol::DrcOptVol()
{
  m_optMat      = 0; 
}
//----------------------------------------------------------------------
DrcOptVol::~DrcOptVol()
{
  delete m_optMat; 
}
//----------------------------------------------------------------------
DrcOptVol* DrcOptVol::clone() const
{
  return new DrcOptVol(*this);
}
//----------------------------------------------------------------------
void DrcOptVol::copy(const DrcOptVol& d)
{
  if (!(&(d.optMaterial())))
    {
      cerr<<" ***  DrcOptVol::copy: optMaterial undefined."<<endl;
      
      exit(EXIT_FAILURE);
    }
  
  if (&(d.optMaterial())) m_optMat = (d.optMaterial()).clone();
}
//----------------------------------------------------------------------
DrcOptVol::DrcOptVol(const DrcOptVol& d) : DrcOptDev(d)
{  
  if (d.m_verbosity>=1) cout<<"  DrcOptVol::DrcOptVol"
			    <<"(const DrcOptVol&) name,copy: "
			    <<d.m_name<<" "<<d.m_copyNumber<<endl;
  copy(d);
}
//----------------------------------------------------------------------
DrcOptVol& DrcOptVol::operator=(const DrcOptVol& d)
{
  if (d.m_verbosity>=1) cout<<"  DrcOptVol::operator="
			    <<"(const DrcOptVol&) name,copy: "
			    <<d.m_name<<" "<<d.m_copyNumber<<endl;
  if (&d != this)
    {
      static_cast<DrcOptDev&>((*this)) = d; // assignment of base class part.
      copy(d);
    }
  return *this;
}
//----------------------------------------------------------------------
void DrcOptVol::setOptMaterial(const DrcOptMatAbs& mat)
{
  if (&mat)
    {
      delete m_optMat;
      m_optMat = mat.clone();
    }
}
//----------------------------------------------------------------------
void DrcOptVol::propagate(DrcPhoton& ph)
{
  static const double eps = 1.0e-9;


  if (verbosity()>=4) cout<<"    DrcOptVol::propagate"<<endl;
  //if (verbosity()>=4) ph.print();
  list<DrcSurfAbs*>::const_iterator isurf;

 
  while (ph.fate() == Drc::PhotFlying)
    {// while
      ph.setDevice(this);
      //---------------
      // search the closest surface which is not the surface where the photon is 
      // right now. 
      DrcSurfAbs* surf_closest=0;
      double path_length_min = numeric_limits<double>::max();
      for(isurf=m_listSurf.begin(); isurf != m_listSurf.end(); ++isurf) 
	{
	  if (verbosity()>=4) cout<<"     check "<<(*isurf)->name()
				  <<" (surface coupling = "
				  <<(*isurf)->coupled()<<")"<<endl;
	  XYZPoint pos_new;
	  double   path_length;
	  if ((*isurf)->surfaceHit(ph,pos_new,path_length))
	    {
	      if (verbosity()>=4) cout<<"     hit for "<<(*isurf)->name()<<endl;
	      if (verbosity()>=4) cout<<"     path_length "<<path_length<<endl;
	      if (path_length<path_length_min && path_length>eps)
		{
		  path_length_min = path_length;
		  surf_closest    = (*isurf);
		}
	    }
	  else
	    {
	      if (verbosity()>=4) cout<<"     no hit "<<(*isurf)->name()<<endl;
	    }

	}
      if (surf_closest==0)
	{
	  if (ph.fate() == Drc::PhotFlying) 
	    {
	      if (verbosity()>=4) 
		cout<<"     no hit with any surface ->edge hit->lost"<<endl;
	      if (m_photonTrace) ph.print(*m_photonTraceStream);
	      ph.setFate(Drc::PhotLost);
	      break;
	    }
	}
      // reset photon, that was only a search...
      ph.setFate(Drc::PhotFlying);

      //---------------


      if (verbosity()>=4) cout<<"     closest "<<surf_closest->name()<<endl;
      XYZPoint pos_new;
      double   path_length;
      if (surf_closest->surfaceHit(ph,pos_new,path_length))
	{// hit surf_closest

	  if (verbosity()>=4) 
	    {
	      XYZPoint pos_old(ph.position());
	      cout<<"     hit, set pos from "<<pos_old<<endl;
	      cout<<"                  to   "<<pos_new<<endl;
	    }

	  //if (ph.fate()!=Drc::PhotFlying) return;

	  // ----------- Absorption 
	  if (optMaterial().absorptionFlag(ph.wavelength(),path_length))
	    {
	      //ph.setPosition(ph.position()+path_length*ph.direction());
	      if (m_photonTrace) ph.print(*m_photonTraceStream);
	      ph.setFate(Drc::PhotAbsorbed);
	      break; // leave while loop
	    }
	  else
	    {
	      ph.setPosition(pos_new);
	      if (m_photonTrace) ph.print(*m_photonTraceStream);
	      if (ph.fate()!=Drc::PhotFlying) break;//###1
	    }
	  //------------

	  //------------ check coupled volumes

	  if (verbosity()>=4) cout<<"     surface coupling = "
				  <<surf_closest->coupled()
				  <<" for "
				  <<name()<<" "
				  <<copyNumber()<<" "
				  <<surf_closest->name()<<endl;

	  if (surf_closest->coupled())
	    {// if surf_closest coupled
	      list<DrcOptDev*>::const_iterator  idev_coupled;
	      list<DrcSurfAbs*>::const_iterator isurf_coupled;
	      
	      for (idev_coupled  = (surf_closest->coupledDeviceList()).begin(),
		     isurf_coupled = (surf_closest->coupledSurfaceList()).begin(); 
		   idev_coupled != (surf_closest->coupledDeviceList()).end(); 
		   ++idev_coupled, ++isurf_coupled)
		{
		  if (verbosity()>=4) cout<<"     coupling to dev,copy,surf= "
					  <<(*idev_coupled)->name()<<","
					  <<(*idev_coupled)->copyNumber()<<","
					  <<(*isurf_coupled)->name()<<endl;

		  // deal with eps problem in sphere and parabolid
		  DrcPhoton ph1(ph);
		  if (! (*isurf_coupled)->isFlat())
		  {
		    ph1.setPosition(ph.position()-ph.direction()*0.1);
		    if (verbosity()>=4) cout<<"     bring back from "
		  		      <<ph.position()<<" to "
		  		      <<ph1.position()<<endl;
		  }

		  bool hit = (*isurf_coupled)->surfaceHit(ph1,pos_new,path_length);

		  if (verbosity()>=4) cout<<"     hit="<<hit<<" with "
					  <<(*isurf_coupled)->name()<<endl;
		  if (hit)
		    {

		      if (&(surf_closest->reflectivity())) // Reflectivity defined
			{
			  // check reflectivity
			  XYZVector norm = surf_closest->normal(ph.position());
			  Drc::Reflectivity refl = Drc::ReflReflected;
			  refl = surf_closest->reflectivity().reflectivity(ph,norm);
		      
		      
			  if (refl == Drc::ReflAbsorbed)
			    {
			      if (verbosity()>=4) 
				cout<<"     DrcOptMirror::propagate: mirror absorbed"<<endl;
			      ph.setFate(Drc::PhotAbsorbed);
			      return;
			    }
			  if (refl == Drc::ReflTransmitted)
			    {
			      // do nothing
			    }
			  if (refl == Drc::ReflReflected)
			    {
			      ph.reflect(norm);
			    }
			  if (refl == Drc::ReflRefracted)
			    {
			      //if (refract(ph,norm))
			      //{
			      //  ph.setFate(Drc::PhotLost); // Photon refracted in nirvana.
			      //  if (verbosity()>=4) cout<<"     Photon lost"<<endl;
			      //}
			      // do nothing
			    }

			}


		      // bring photon inside coupled volume
		      // to prevent infinite recursion of coupled surface hits.
		      // The factor 2 comes from comparisons with same eps in 
		      // surfaceHit.

		      DrcOptMatAbs* opt_mat = &((*idev_coupled)->optMaterial());
		      if (verbosity()>=4) cout<<"     opt_mat="<<opt_mat<<endl;
		      if (opt_mat)
			{
			  double n2 = opt_mat->refIndex(ph.wavelength());
			  bool iref = refract(ph,
					      surf_closest->normal(ph.position()),
					      n2);

			  if (verbosity()>=4) cout<<" refract in new volume flag = "
						  <<iref<<endl;

			  if (iref) // refraction in new volume
			    {
			      if (verbosity()>=4) cout<<"     go into new volume "
						      <<(*idev_coupled)->name()<<endl;
			      (*idev_coupled)->propagate(ph);
			    }
			  else
			    {
			      break; // coupled surface loop, since phot. is reflected.
			    }
			}
		      else // screen or mirror...
			{
			  // do not bring photons inside flat objects
			  // ph.setPosition(ph.position()+2*eps*ph.direction());
			  (*idev_coupled)->propagate(ph);
			}
		      if (ph.fate() != Drc::PhotFlying) break;
		    }
		}
	    }
	  //----------
	  else
	    {// uncoupled surfaces
	      XYZVector norm = surf_closest->normal(ph.position());

	      if (&(surf_closest->reflectivity())) // Reflectivity defined
		{
		  ph.setDevice(this);

		  Drc::Reflectivity refl = Drc::ReflReflected;
		  refl = surf_closest->reflectivity().reflectivity(ph,norm);
			   

		  if (refl == Drc::ReflAbsorbed)
		    {
		      if (verbosity()>=4) 
			cout<<"     DrcOptMirror::propagate: mirror absorbed"<<endl;
		      ph.setFate(Drc::PhotAbsorbed);
		      return;
		    }
		  if (refl == Drc::ReflTransmitted)
		    {
		      // do nothing
		    }
		  if (refl == Drc::ReflReflected)
		    {
		      ph.reflect(norm);
		    }
		  if (refl == Drc::ReflRefracted)
		    {
		      if (refract(ph,norm))
			{
			  ph.setFate(Drc::PhotLost); // Photon refracted in nirvana.
			  if (verbosity()>=4) cout<<"     Photon lost"<<endl;
			}
		    }
		}
	      else if (refract(ph,norm) && ph.fate() == Drc::PhotFlying) 
		{//perfect surface
		  ph.setFate(Drc::PhotLost); // Photon refracted in nirvana.
		  if (verbosity()>=4) cout<<"     Photon lost"<<endl;
		}
	      // else, do nothing, photon is reflected.
	    }
	  if (ph.fate() != Drc::PhotFlying) break; // while loop
	}// if surf_closest... 
    }
}
//----------------------------------------------------------------------
bool DrcOptVol::refract( DrcPhoton& ph, XYZVector normal, double n2) const
{
  static const double eps = 1.0e-9;

  bool     refract_flag;

  XYZVector dir1   = ph.direction();

  double   lambda = ph.wavelength();
  double   n1     = optMaterial().refIndex(lambda);

  // norm points outside.
  XYZVector norm = (dir1.Dot(normal)>0) ? normal : -normal;



  //
  // sin(alpha1)   n2
  // ----------- = --  (Snellius)
  // sin(alpha2)   n1
  //

  // alpaha1, alpha2 > 0

  if (verbosity()>=4) cout<<"     norm="
			  <<norm.X()<<" "<<norm.Y()<<" "<<norm.Z()<<endl;
  if (verbosity()>=4) cout<<"     dir1="
			  <<dir1.X()<<" "<<dir1.Y()<<" "<<dir1.Z()<<endl;
  if (verbosity()>=4) cout<<"     norm*dir="<<norm.Dot(dir1)<<endl;

  double alpha1 = acos(norm.Dot(dir1));
  double alpha2;  
  if (sin(alpha1) * n1/n2 > 1) // reflect photon
    {
      refract_flag = false;
      ph.reflect(norm);
      //ph.setDirection(dir1 - 2*(norm.Dot(dir1)*norm)); // new direction
      if (verbosity()>=4) cout<<"     reflection"<<endl;
    }
  else                         // refract photon
    {
      refract_flag = true;

      if (verbosity()>=4) cout<<"     alpha1,n1,n2 = "<<alpha1<<" "<<n1<<" "<<n2<<endl;
      alpha2 = asin(sin(alpha1) * n1/n2);
    
      XYZVector rhelp = norm.Cross(dir1);
      XYZVector diff  = rhelp.Cross(norm);
      diff = diff.Unit();
      
      // diff perp to norm points in direction to the tip of vel.
      //
      // dir2 = norm + x*diff such that L(x*diff)/L(n) = tan(alpha2)
      //
      double x = tan(alpha2);
      if (verbosity()>=4) cout<<"     x,alpha2 = "<<x<<" "<<alpha2<<endl;
      XYZVector dir2 = norm + x*diff;

      if (verbosity()>=4) cout<<"     dir2 before unit = "<<
	dir2.X()<<" "<<
	dir2.Y()<<" "<<
	dir2.Z()<<endl; 
      dir2 = dir2.Unit();
      ph.setDirection(dir2);
      ph.setPosition(ph.position() + 2*eps*dir2); // bring it in new volume or outside.
      if (verbosity()>=4) 
	{
	  cout<<"     refraction, brought it inside to "
	      <<ph.position().X()<<" "
	      <<ph.position().Y()<<" "
	      <<ph.position().Z()<<endl;
	  cout<<"     refraction, brought in direction "
	      <<ph.direction().X()<<" "
	      <<ph.direction().Y()<<" "
	      <<ph.direction().Z()<<endl;
	}
    } 

  return refract_flag;

}