#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 <limits>

#include <fstream>
using std::fstream;

#include <utility>
using std::pair;

#include <map>
using std::map;

//#include <cmath>

#include "TVector3.h"
#include "TRandom.h"
#include "TRotation.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 "PndDrcSurfPolyAsphere.h"
#include "PndDrcSurfPolySphere.h"
#include "PndDrcSurfPolyPara.h"
#include "PndDrcSurfPolyAsphere.h"
#include "PndDrcOptReflSilver.h"
#include "PndDrcOptReflNone.h"
#include "PndDrcOptReflPerfect.h"
#include "PndDrcOptMatAbs.h"
#include "PndDrcOptMatLithotecQ0.h"
#include "PndDrcOptMatNLAK33A.h"
#include "PndDrcOptMatLLF1.h"
#include "PndDrcOptMatVacuum.h"
#include "PndDrcOptDev.h"
#include "PndDrcOptDevSys.h"
#include "PndDrcOptVol.h"
#include "PndDrcOptDevManager.h"
#include "PndDrcUtil.h"


int main(int argc, char *argv[])
{
  //double conical_const =  0; //sphere
  double conical_const = 0.5; //parabola

  //int ioption = 1; // parallel beam (0 deg)
  //int ioption = 2; // 0,  +- 20,  +-40 deg
  //int ioption = 3; // C-cone
  //int ioption = 4; // random in front of lens.
  //int ioption = 5; // single photon for debugging.
  int ioption = 6; // grid 5 deg

  double pi=3.1415926535;

  PndDrcOptReflPerfect refl_perfect;
  PndDrcOptReflNone    refl_none;
  double            lambda=500;

  double slab_width  = 17;
  double slab_height = 17;


  PndDrcOptMatLithotecQ0 quartz;
  PndDrcOptMatVacuum     vacuum;


  // create a lens with slab lateral dimensions

  XYZPoint q0(-slab_width/2, +slab_height/2, 0.0);
  XYZPoint q1(-slab_width/2, -slab_height/2, 0.0);
  XYZPoint q2(+slab_width/2, -slab_height/2, 0.0);
  XYZPoint q3(+slab_width/2, +slab_height/2, 0.0);

  PndDrcSurfPolyAsphere lens_sphere;
  lens_sphere.setVerbosity(0);
  lens_sphere.addPoint(q0);
  lens_sphere.addPoint(q1);
  lens_sphere.addPoint(q2);
  lens_sphere.addPoint(q3);
  lens_sphere.setRadius(95); // f = R/(n-1) = 200 mm
  lens_sphere.setPrintColor(3);
  lens_sphere.setConicalConstant(conical_const);
  lens_sphere.setName("lens_sphere");
  double     dist = lens_sphere.centerPoint().Z();
  lens_sphere.addTransform(Transform3D(XYZVector(0,0,-(dist))));

  q0+=XYZVector(0,0,1.5);
  q1+=XYZVector(0,0,1.5);
  q2+=XYZVector(0,0,1.5);
  q3+=XYZVector(0,0,1.5);

  PndDrcSurfPolyFlat lens_base;
  lens_base.setVerbosity(0);
  lens_base.addPoint(q0);
  lens_base.addPoint(q1);
  lens_base.addPoint(q2);
  lens_base.addPoint(q3);
  lens_base.setPrintColor(3);
  lens_base.setName("lens_base");

  PndDrcSurfQuadFlatDiff lens_side1;
  lens_side1.addSurface(lens_base,  q0,q1);
  lens_side1.addSurface(lens_sphere,
			lens_sphere.limitingPoint(0),lens_sphere.limitingPoint(1));
  lens_side1.setPrintColor(3);
  lens_side1.setName("lens_side1");
  PndDrcSurfQuadFlatDiff lens_side2;
  lens_side2.setVerbosity(0);
  lens_side2.addSurface(lens_base,  q1,q2);
  lens_side2.addSurface(lens_sphere,
			lens_sphere.limitingPoint(1),lens_sphere.limitingPoint(2));
  lens_side2.setPrintColor(3);
  lens_side2.setName("lens_side2");
  PndDrcSurfQuadFlatDiff lens_side3;
  lens_side3.addSurface(lens_base,  q2,q3);
  lens_side3.addSurface(lens_sphere,
			lens_sphere.limitingPoint(2),lens_sphere.limitingPoint(3));
  lens_side3.setPrintColor(3);
  lens_side3.setName("lens_side3");
  PndDrcSurfQuadFlatDiff lens_side4;
  lens_side4.addSurface(lens_base,  q3,q0);
  lens_side4.addSurface(lens_sphere,
			lens_sphere.limitingPoint(3),lens_sphere.limitingPoint(0));
  lens_side4.setPrintColor(3);
  lens_side4.setName("lens_side4");


  lens_side1.setReflectivity(refl_none);
  lens_side2.setReflectivity(refl_none);
  lens_side3.setReflectivity(refl_none);
  lens_side4.setReflectivity(refl_none);
  

  PndDrcOptVol lens;
  lens.setVerbosity(0);
  lens.addSurface(lens_base);
  lens.addSurface(lens_sphere);
  lens.addSurface(lens_side1);
  lens.addSurface(lens_side2);
  lens.addSurface(lens_side3);
  lens.addSurface(lens_side4);
  lens.setOptMaterial(vacuum);
  lens.setName("lens");
  lens.addTransform(Transform3D(XYZVector(0,0,2)));




  XYZPoint p1(-300.0, +300.0, 0);  //            p5----------p8
  XYZPoint p2(-300.0, -300.0, 0);  //           /|          /|
  XYZPoint p3(+300.0, -300.0, 0);  //          / |         / |
  XYZPoint p4(+300.0, +300.0, 0);  //         /  p6-------/--p7
  //                                              /  /        /  /
  //                                             /  /        /  /
  //                                            /  /        /  /
  //                                           /  /        /  /
  XYZPoint p5(-300.0, +300.0, +300); // p1-s1--s2---p4 /
  XYZPoint p6(-300.0, -300.0, +300); // | /s5--s6   | /
  XYZPoint p7(+300.0, -300.0, +300); // |/ s7--s8   |/
  XYZPoint p8(+300.0, +300.0, +300); // p2-s3--s4--p3 

  PndDrcSurfPolyFlat box_side1;
  box_side1.addPoint(p1);
  box_side1.addPoint(p2);
  box_side1.addPoint(p3);
  box_side1.addPoint(p4);
  box_side1.setReflectivity(refl_none);
  box_side1.setName("box_side1");
  PndDrcSurfPolyFlat box_side2;
  box_side2.addPoint(p2);
  box_side2.addPoint(p3);
  box_side2.addPoint(p7);
  box_side2.addPoint(p6);
  box_side2.setReflectivity(refl_none);
  box_side2.setName("box_side2");
  PndDrcSurfPolyFlat box_side3;
  box_side3.addPoint(p3);
  box_side3.addPoint(p4);
  box_side3.addPoint(p8);
  box_side3.addPoint(p7);
  box_side3.setReflectivity(refl_none);
  box_side3.setName("box_side3");
  PndDrcSurfPolyFlat box_side4;
  box_side4.addPoint(p1);
  box_side4.addPoint(p4);
  box_side4.addPoint(p8);
  box_side4.addPoint(p5);
  box_side4.setReflectivity(refl_none);
  box_side4.setName("box_side4");
  PndDrcSurfPolyFlat box_side5;
  box_side5.addPoint(p5);
  box_side5.addPoint(p6);
  box_side5.addPoint(p7);
  box_side5.addPoint(p8);
  box_side5.setPixel();
  box_side5.setName("box_side5");
  PndDrcSurfPolyFlat box_side6;
  box_side6.addPoint(p5);
  box_side6.addPoint(p6);
  box_side6.addPoint(p2);
  box_side6.addPoint(p1);
  box_side6.setPixel();
  box_side6.setName("box_side6");



  PndDrcOptVol box;
  box.setVerbosity(0);
  box.addSurface(box_side1);
  box.addSurface(box_side2);
  box.addSurface(box_side3);
  box.addSurface(box_side4);
  box.addSurface(box_side5);
  box.addSurface(box_side6);
  box.setOptMaterial(quartz);
  box.setName("box");

  XYZPoint s5(-slab_width/2, +slab_height/2, 0);
  XYZPoint s6(+slab_width/2, +slab_height/2, 0);
  XYZPoint s7(-slab_width/2, -slab_height/2, 0);
  XYZPoint s8(+slab_width/2, -slab_height/2, 0);
  XYZPoint t5(-slab_width/2, +slab_height/2, -1000);
  XYZPoint t6(+slab_width/2, +slab_height/2, -1000);
  XYZPoint t7(-slab_width/2, -slab_height/2, -1000);
  XYZPoint t8(+slab_width/2, -slab_height/2, -1000);
  //s 5 6 at dist-1
  //  7 8
  //
  //t 5 6 at dist-1-1000
  //  7 8
  //
  PndDrcSurfPolyFlat slab_side1;
  slab_side1.setVerbosity(0);
  slab_side1.addPoint(s5); // face
  slab_side1.addPoint(s6);
  slab_side1.addPoint(s8);
  slab_side1.addPoint(s7);
  slab_side1.setName("slab_side1");
  PndDrcSurfPolyFlat slab_side2;
  slab_side2.addPoint(s5);//left
  slab_side2.addPoint(s7);
  slab_side2.addPoint(t7);
  slab_side2.addPoint(t5);
  slab_side2.setName("slab_side2");
  PndDrcSurfPolyFlat slab_side3;
  slab_side3.addPoint(s7);//down
  slab_side3.addPoint(s8);
  slab_side3.addPoint(t8);
  slab_side3.addPoint(t7);
  slab_side3.setName("slab_side3");
  PndDrcSurfPolyFlat slab_side4;
  slab_side4.addPoint(s6);//right
  slab_side4.addPoint(s8);
  slab_side4.addPoint(t8);
  slab_side4.addPoint(t6);
  slab_side4.setName("slab_side4");
  PndDrcSurfPolyFlat slab_side5;
  slab_side5.addPoint(s5);//top
  slab_side5.addPoint(s6);
  slab_side5.addPoint(t6);
  slab_side5.addPoint(t5);
  slab_side5.setName("slab_side5");
  PndDrcSurfPolyFlat slab_side6;
  slab_side6.addPoint(t5); // face
  slab_side6.addPoint(t6);
  slab_side6.addPoint(t8);
  slab_side6.addPoint(t7);
  slab_side6.setName("slab_side6");

  //slab_side1.setReflectivity(refl_perfect);
  slab_side2.setReflectivity(refl_perfect);
  slab_side3.setReflectivity(refl_perfect);
  slab_side4.setReflectivity(refl_perfect);
  slab_side5.setReflectivity(refl_perfect);
  slab_side6.setReflectivity(refl_perfect);

  PndDrcOptVol slab;
  slab.setVerbosity(0);
  slab.addSurface(slab_side1);
  slab.addSurface(slab_side2);
  slab.addSurface(slab_side3);
  slab.addSurface(slab_side4);
  slab.addSurface(slab_side5);
  slab.addSurface(slab_side6);
  slab.setOptMaterial(quartz);
  slab.setName("slab");




  PndDrcOptDevSys opt_system;
  opt_system.setName("opt_system");
  opt_system.setVerbosity(0);
  opt_system.addDevice(box);
  opt_system.addDevice(slab);
  opt_system.embedDevice(lens,box);

  // couple surface 1 of device 1 with surface 2 of device 2
  //                       dev1   dev2    surf1     surf2
  opt_system.coupleDevice("slab","box","slab_side1", "box_side1");




  PndDrcOptDevManager* manager = new PndDrcOptDevManager();
  manager->addDeviceSystem(opt_system);


  fstream geo;
  geo.open("Geo.C",std::ios::out);
  geo<<"{"<<endl;
  geo<<"    TCanvas *c1 = new TCanvas(\"c1\"); "<<endl;
  geo<<"    TView *view = new TView(1);"<<endl;
  //geo<<"    view->SetRange(-450,-450,-300,450,450,500);"<<endl;
   geo<<"    view->SetRange(-450,-450,-900,450,450,900);"<<endl;
 geo<<"    Int_t i;"<<endl;
  geo<<"    view->SetView(180,45,90,i);"<<endl;
  geo<<"    view->SetView(90,0,90,i);"<<endl;
  geo<<"    view->Zoom();"<<endl;
  geo<<"    view->Zoom();"<<endl;
  geo<<"    view->Zoom();"<<endl;
  geo<<"    view->Zoom();"<<endl;
  geo<<"    view->Zoom();"<<endl;
  geo<<"    view->Zoom();"<<endl;
  geo<<"    view->Zoom();"<<endl;
  
  // the following command sets a flag within the manager and all photons from
  // now on will be traced and can be plotted by calling within root
  // .x Geo.C 
  // .x Screen.C
  // 
  manager->print(geo);
  

  //geo<<"}"<<endl;
  //exit(0);

  list<PndDrcPhoton> list_photon; // get list

  PndDrcPhoton ph;


  bool photons_exist = false;
  if (ioption==1)  
    {
      slab_width*=0.5;
      slab_height*=0.5;

      for (double xx=-slab_width/2; xx<=slab_width/2; xx+=slab_width/4)
	{
	  for (double yy=-slab_height/2; yy<=slab_height/2; yy+=slab_height/4)
	    {
	      ph.setPosition(XYZPoint(xx,yy,dist-500));
	      ph.setDirection(XYZVector(0,0,1));
	      ph.setWavelength(lambda);
	      list_photon.push_back(ph);
	      //goto raus;	      
	    }
	}
      // raus:
      photons_exist = true;
      manager->setPhotonList(list_photon,"slab","opt_system",0,0);
    }
  if (ioption==2)  
    {
      for (double angle1=-30; angle1<30.5; angle1+=5.0)
	{
	  for (double xx=-0.2*slab_width; xx<=0.21*slab_width; xx+=.40*slab_width)
	    {
	      double yy=0;
	      //for (double lambda1=330; lambda1<631; lambda1+=100)
	      {
		double angle;
		double lambda1 = lambda;
		ph.setPosition(XYZPoint(xx,yy,-500.5));
	      
		angle = pi/180.0 * angle1;
		ph.setDirection(XYZVector(sin(+angle),0,cos(+angle)));
		ph.setWavelength(lambda1);
		list_photon.push_back(ph);
	     
		   
	      }
	    }
	}
      photons_exist = true;
      manager->setPhotonList(list_photon,"slab");
    }
  if (ioption==3)
    {
      XYZPoint pos(0,-10,-45.0);
      XYZVector dir(0.0,1,-3); 
      double   beta = 0.80;
      photons_exist = manager->cerenkov(pos,dir,beta,100,200); // generate photons
    }
  if (ioption==4)  
    {
      TRandom ran;
      for (int ii=0; ii<30; ii++)
	{
	  ph.setPosition(XYZPoint(0,0,dist-20));
	  double phi =   ran.Uniform(0,2*pi);
	  double costh = ran.Uniform(-1,1);
	  double sinth = sqrt(1-costh*costh);
	  ph.setDirection(XYZVector(sinth*sin(phi),sinth*cos(phi),costh));
	  ph.setWavelength(lambda);
	  list_photon.push_back(ph);
	}
	
      photons_exist = true;
      manager->setPhotonList(list_photon,"slab");
    }
  if (ioption==5)
    {
      double angle = pi/180.0 * 0;
      ph.setDirection(XYZVector(0,sin(+angle),cos(+angle)));

      ph.setPosition(XYZPoint(0,5,-5));
      ph.setWavelength(400);
      list_photon.push_back(ph);
      ph.setPosition(XYZPoint(0,-5,-5));
      ph.setWavelength(400);
      list_photon.push_back(ph);


      ph.setPosition(XYZPoint(0,5,-5));
      ph.setWavelength(600);
      list_photon.push_back(ph);
      ph.setPosition(XYZPoint(0,-5,-5));
      ph.setWavelength(600);
      list_photon.push_back(ph);



      photons_exist = true;
      manager->setPhotonList(list_photon,"slab");
    }
  if (ioption==6)  
    {
      // straight lines.
      TRandom ran;
      for (double angle=5; angle<=10.5; angle+=5)
	{

	  cout<<" angle = "<<angle<<endl;
	  for (double lambda1=630; lambda1>329; lambda1-=10)
	    {
	      //double lambda1=lambda;
	      // go in x dir
	      double y1 = tan(angle*pi/180);
	      double x1 = y1;
	      double scale = 3/angle;
	      for (double y=-y1; y<= y1; y+= 2*y1/50*scale)
		{
		  //double xx=ran.Uniform(-0.5*slab_width,0.5*slab_width);
		  //double yy=ran.Uniform(-0.5*slab_height,0.5*slab_height);
		  ph.setPosition(XYZPoint(0,0,-0.5));
		  double x = x1;
		  ph.setDirection(XYZVector(x,y,1).Unit());
		  ph.setWavelength(lambda1);
		  list_photon.push_back(ph);
		  ph.setDirection(XYZVector(-x,y,1).Unit());
		  ph.setWavelength(lambda1);
		  list_photon.push_back(ph);
		}
	      for (double x=-x1; x<= x1; x+= 2*y1/50*scale)
		{
		  //double xx=ran.Uniform(-0.5*slab_width,0.5*slab_width);
		  //double yy=ran.Uniform(-0.5*slab_height,0.5*slab_height);
		  ph.setPosition(XYZPoint(0,0,-0.5));
		  double y = y1;
		  ph.setDirection(XYZVector(x,y,1).Unit());
		  ph.setWavelength(lambda1);
		  list_photon.push_back(ph);
		  ph.setDirection(XYZVector(x,-y,1).Unit());
		  ph.setWavelength(lambda1);
		  list_photon.push_back(ph);
		}
	     
	    }
	}
      photons_exist = true;
      manager->setPhotonList(list_photon,"slab");
    }

  if (photons_exist)
    {
      manager->propagate();                           // propagate photons
    }

  list_photon = manager->photonList(); // get list




  geo<<"}"<<endl;
  //
  // the intention is to play around with routines.

  // create a list of photons in bar
  //exit(1);
  geo.close();
  
  // write to screen...

  fstream scr;
  scr.open("Screen.C",std::ios::out);
  scr<<"{"<<endl;
  scr<<"    TCanvas *c1 = new TCanvas(\"c1\"); "<<endl;
  scr<<"    TH1F *hgr = new TH1F(\"hgr\",\"spatial position [mm]\",600,-300,300);"<<endl;
  scr<<"    hgr->SetMarkerStyle(7);"<<endl;
  scr<<"    hgr->SetMarkerSize(0.5);"<<endl;
  scr<<"    hgr->SetMinimum(-300);"<<endl;
  scr<<"    hgr->SetMaximum(+300);"<<endl;
  scr<<"    hgr->Draw(\"POL\");"<<endl;
  int icnt_measured = 0;
  int icnt_flying   = 0;
  int icnt_lost     = 0;
  int icnt_absorbed = 0;

  //TRotation rotInv_screen = rot_screen.Inverse();
  list<PndDrcPhoton>::iterator iph;
  PndDrcPhoton ph_old;

  double dista,distb;
  for(iph=list_photon.begin(); iph != list_photon.end(); ++iph) 
    {


      //(*iph).print();
      if ((*iph).fate()==Drc::PhotMeasured)
	{
	  icnt_measured++;
	  XYZPoint pos = (*iph).position();
	  //pos -= shift_screen;
	  //pos.Transform(rotInv_screen);

	  double x     = pos.X();
	  double y     = pos.Y();

	  //cout<<" i,p,d="<<icnt_measured<<" "
	  //<<(*iph).position()<<" "<<(*iph).direction()<<endl;

	  if ((icnt_measured%2)==0)
	    {
	      //cout<<icnt_measured<<endl;
	      XYZPoint p1((*iph).position());
	      XYZVector d1((*iph).direction());
	      XYZPoint p2(ph_old.position());
	      XYZVector d2(ph_old.direction());
	      double d1x = d1.X();
	      double d1z = d1.Z();
	      double d2x = d2.X();
	      double d2z = d2.Z();
	      double p1x = p1.X();
	      double p1z = p1.Z();
	      double p2x = p2.X();
	      //double p2z = p2.Z();


	      double d = (d1x - d1z/d2z*d2x);
	      double lam = (p2x-p1x)/d;
	      //cout<<p1x<<" "<<p1z<<" "<<d1x<<" "<<d1z<<endl;
	      //cout<<p2x<<" "<<p2z<<" "<<d2x<<" "<<d2z<<endl;
	      
	      //cout<<" "<<(p2x-p1x)<<" "<<d<<" "<<lam*d1z<<endl;
	      //cout<<" distz = "<<p1z+lam*d1z<<endl;
	      //cout<<" distx = "<<p1x+lam*d1x<<endl;

	      if (icnt_measured==2) dista = p1z+lam*d1z;
	      if (icnt_measured==14) distb = p1z+lam*d1z;

	    }

	  scr<<"    TMarker* t = new TMarker("<<x<<","<<y<<",20);"<<endl;
	  scr<<"    t->SetMarkerStyle(7);"<<endl;
	  scr<<"    t->SetMarkerColor("
		<<(*iph).colorNumber((*iph).wavelength())
		<<");"<<endl;
	  scr<<"    t->SetMarkerSize(0.7);"<<endl;
	  scr<<"    t->Draw();"<<endl;
	  ph_old = (*iph);
	}
      else if ((*iph).fate()==Drc::PhotFlying)   icnt_flying++; // should never happen.
      else if ((*iph).fate()==Drc::PhotAbsorbed) icnt_absorbed++;
      else                                       
	{
	  icnt_lost++;
	  //XYZPoint pos = (*iph).position();
	  //cout<<pos<<endl;
	}
    }
  //cout<<" ratio = "<<dista/distb<<endl;
  scr<<"}"<<endl;
  scr.close();

  int icnt = icnt_measured+icnt_flying+icnt_lost+icnt_absorbed;
  cout<<" generated photons: "<<icnt<<endl;
  cout<<" measured  photons: "<<icnt_measured<<endl;
  cout<<" absorbed  photons: "<<icnt_absorbed<<endl;
  cout<<" lost      photons: "<<icnt_lost<<endl;
  delete manager;

  return EXIT_SUCCESS;

}