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

#include <valarray>
using std::valarray;

#include <vector>
using std::vector;

#include <string>
using std::string;

#include <list>
using std::list;


#include <map>
using std::map;

#include <limits>

#include <fstream>
using std::fstream;

#include <utility>
using std::pair;
//#include <cmath>

#include "TObject.h"
#include "TVector3.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 "PndDrcOptMatAbs.h"
#include "PndDrcOptDev.h"
#include "PndDrcUtil.h"


//----------------------------------------------------------------------
PndDrcPhoton::PndDrcPhoton()
{
  m_lambda      = 0;
  m_position    = XYZPoint(0,0,0);
  m_positionOld = XYZPoint(0,0,0);
  m_direction   = XYZPoint(0,0,0);
  m_fate        = Drc::PhotFlying;
  m_reflections = 0;
  m_verbosity   = 0;
  m_time        = 0;
  m_dev         = 0;
}
//----------------------------------------------------------------------
void PndDrcPhoton::setPosition(const XYZPoint& pos)
{
  if (m_dev && m_dev->radiator()) // no flat device
    {
      double n = (m_dev->optMaterial()).refIndex(m_lambda);
      double dndl = (m_dev->optMaterial()).refIndexDeriv(m_lambda);
      double len = sqrt((pos-m_position).Mag2()); // mm
      double v_phase = 299.792/n;
      double v_group =  v_phase * ( 1.0 - m_lambda/n*dndl);
      double time = len/v_group;
      m_time += time;          // pos in mm time in ns
    }

  m_positionOld = m_position;
  m_position    = pos;
};
//----------------------------------------------------------------------
bool PndDrcPhoton::refract(XYZVector normal, double n1, double n2) 
{
  static const double eps = 1.0e-9;

  bool      refract_flag;
  XYZVector dir1   = direction();
 
 // 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;
      cout<<"     dir1="<<dir1.X()<<" "<<dir1.Y()<<" "<<dir1.Z()<<endl;
      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;
      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();
      setDirection(dir2);
      setPosition(position() + 2*eps*dir2); // bring it in new volume or outside.
      if (verbosity()>=4) 
	{
	  cout<<"     refraction, brought it inside to "
	      <<position().X()<<" "
	      <<position().Y()<<" "
	      <<position().Z()<<endl;
	  cout<<"     refraction, brought in direction "
	      <<direction().X()<<" "
	      <<direction().Y()<<" "
	      <<direction().Z()<<endl;
	}
    } 

  return refract_flag;

}

//----------------------------------------------------------------------
void PndDrcPhoton::reflect(const XYZVector& normal)
{
  XYZVector normal1(normal);

  if (m_direction.Dot(normal) < 0) normal1 *= -1;
  m_direction = m_direction-2*(m_direction.Dot(normal1))*normal1;

  m_reflections++;

}
//----------------------------------------------------------------------
int PndDrcPhoton::colorNumber(double lambda) const
{

  // the colors were taken from wikipedia entry about colors.

  if      (lambda< 0.1) return 17; // grey
  if      (lambda> 740) return 50; // dark red
  else if (lambda> 625) return 2;  // red
  else if (lambda> 590) return 42; // orange
  else if (lambda> 565) return 5;  // yellow
  else if (lambda> 520) return 3;  // green
  else if (lambda> 500) return 7;  // cyan
  else if (lambda> 450) return 4;  // blue
  else if (lambda> 430) return 9;  // indigo
  else if (lambda> 380) return 39; // violet

  return 33; // blue grey
}
//----------------------------------------------------------------------
void PndDrcPhoton::print(fstream& stream) const
{

  //if (m_reflections==0) return;

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

  stream<<" l->SetPoint("<<0<<","
	<<m_positionOld.X()<<","
	<<m_positionOld.Y()<<","
	<<m_positionOld.Z()<<");"<<endl;
  stream<<" l->SetPoint("<<1<<","
	<<m_position.X()<<","
	<<m_position.Y()<<","
	<<m_position.Z()<<");"<<endl;
  stream<<" l->SetLineColor("<<colorNumber(m_lambda)<<");"<<endl;
  stream<<" l->Draw();"<<endl;

}
//----------------------------------------------------------------------
void PndDrcPhoton::print() const
{
  cout<<" --- photon ---"<<endl;
  cout<<"     pos    = "<<m_position<<endl;
  cout<<"     dir    = "<<m_direction<<endl;
  cout<<"     time   = "<<m_time<<endl;
  cout<<"     lambda = "<<m_lambda<<endl;
  cout<<"     refl   = "<<m_reflections<<endl;
  string str="unknown";
  if (m_fate==Drc::PhotFlying)   str="flying";
  if (m_fate==Drc::PhotMeasured) str="measured";
  if (m_fate==Drc::PhotLost)     str="lost";
  if (m_fate==Drc::PhotAbsorbed) str="absorbed";
  cout<<"     fate   = "<<str<<endl;
  if (m_dev)
    {
      cout<<"     dev    = "<<m_dev->name()<<endl;
    }
  else
    {
      cout<<"     dev    = "<<"zero pointer"<<endl;
    }
  cout<<" --------------"<<endl;
}