#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 "PndDrcSurfPolySphere.h"
#include "PndDrcSurfPolyPara.h"
#include "PndDrcOptReflSilver.h"
#include "PndDrcOptReflNone.h"
#include "PndDrcOptReflPerfect.h"
#include "PndDrcOptMatAbs.h"
#include "PndDrcOptMatLithotecQ0.h"
#include "PndDrcOptDev.h"
#include "PndDrcOptDevSys.h"
#include "PndDrcOptVol.h"
#include "PndDrcOptDevManager.h"
#include "PndDrcUtil.h"

int main()
{
  double pi=3.1415926535;
  //int ioption = 1; // parallel beam (0 deg)
  //int ioption = 2; // 0,  +- 20,  +-40 deg
  int ioption = 1; // C-cone

  //int imirror=1; // spherical
  int imirror=2; // paraboloid


  // Example for a simple bar with focussing mirror.

  int verbosity = 0; // 0=quiet, 1=constructors,2=member,3=functionality
  //                    4=photons, 5=everything.

  PndDrcOptReflPerfect refl_perfect;
  PndDrcOptReflNone    refl_none;

  // create a 10x10x10 cm3 cube to serve as scale
  XYZPoint s1(-50.0,150.0, -600);  //         p5----------p8
  XYZPoint s2(-50.0,+50.0, -600);  //        /|          /|
  XYZPoint s3(+50.0,+50.0, -600);  //       / |         / |
  XYZPoint s4(+50.0,150.0, -600);  //      /  p6-------/--p7
  XYZPoint s5(-50.0,150.0, -700);  // p1---------p4  /
  XYZPoint s6(-50.0,+50.0, -700);  // | /         | /
  XYZPoint s7(+50.0,+50.0, -700);  // |/          |/
  XYZPoint s8(+50.0,150.0, -700);  // p2---------p3 

  PndDrcSurfPolyFlat c1,c2,c3,c4,c5,c6;
  c1.addPoint(s1);
  c1.addPoint(s2);
  c1.addPoint(s3);
  c1.addPoint(s4);
  c1.setName("cdown");
  c2.addPoint(s2);
  c2.addPoint(s6);
  c2.addPoint(s7);
  c2.addPoint(s3);
  c2.setName("cside1");
  c3.addPoint(s1);
  c3.addPoint(s5);
  c3.addPoint(s6);
  c3.addPoint(s2);
  c3.setName("cside2");
  c4.addPoint(s4);
  c4.addPoint(s3);
  c4.addPoint(s7);
  c4.addPoint(s8);
  c4.setName("cside3");
  c5.addPoint(s1);
  c5.addPoint(s4);
  c5.addPoint(s8);
  c5.addPoint(s5);
  c5.setName("cside4");
  c6.addPoint(s8);
  c6.addPoint(s7);
  c6.addPoint(s6);
  c6.addPoint(s5);
  c6.setName("cup");
  PndDrcOptVol cube;
  cube.addSurface(c6);
  cube.addSurface(c2);
  cube.addSurface(c3);
  cube.addSurface(c4);
  cube.addSurface(c5);
  cube.addSurface(c1); // colored objects last to prevent black lines on top
  cube.setName("cube");




  // 8 points define a bar
  
  XYZPoint p1(-10.0, +5.0, 1100);  //         p5----------p8
  XYZPoint p2(-10.0, -5.0, 1100);  //        /|          /|
  XYZPoint p3(+10.0, -5.0, 1100);  //       / |         / |
  XYZPoint p4(+10.0, +5.0, 1100);  //      /  p6-------/--p7
  //                                      /  /        /  /
  //                                     /  /        /  /
  //                                    /  /        /  /
  //                                   /  /        /  /
  XYZPoint p5(-10.0, +5.0, -800);  // p1---------p4  /
  XYZPoint p6(-10.0, -5.0, -800);  // | /         | /
  XYZPoint p7(+10.0, -5.0, -800);  // |/          |/
  XYZPoint p8(+10.0, +5.0, -800);  // p2---------p3 

  // Define from points 6 surfaces of the bar. The points have to be given in 
  // the sequence going around the surface, clock- or counterclock-wise.
  // There are 2 additional surfaces, a mirror and a screen.

  // How to produce surfaces by shift and rotate operation is for sake of clearness
  // not shown here, but in one of the other examples.

  // Declare flat surfaces with arbitrary number of points.
  PndDrcSurfPolyFlat a1,a2,a3,a4,a5,a6;

  a1.setVerbosity(verbosity);
  a1.addPoint(p1);
  a1.addPoint(p2);
  a1.addPoint(p3);
  a1.addPoint(p4);
  a1.setName("adown");
  a1.setReflectivity(refl_perfect);
  a1.setPrintColor(4);

  a2.setVerbosity(verbosity);
  a2.addPoint(p2);
  a2.addPoint(p6);
  a2.addPoint(p7);
  a2.addPoint(p3);
  a2.setName("aside1");

  a3.setVerbosity(verbosity);
  a3.addPoint(p1);
  a3.addPoint(p5);
  a3.addPoint(p6);
  a3.addPoint(p2);
  a3.setName("aside2");

  a4.setVerbosity(verbosity);
  a4.addPoint(p4);
  a4.addPoint(p3);
  a4.addPoint(p7);
  a4.addPoint(p8);
  a4.setName("aside3");

  a5.setVerbosity(verbosity);
  a5.addPoint(p1);
  a5.addPoint(p4);
  a5.addPoint(p8);
  a5.addPoint(p5);
  a5.setName("aside4");

  a6.setVerbosity(verbosity);
  a6.addPoint(p8);
  a6.addPoint(p7);
  a6.addPoint(p6);
  a6.addPoint(p5);
  a6.setName("aup");

  // create a volume consiting of surfaces
  // create a material the bar will consist of
  PndDrcOptVol bar;
  PndDrcOptMatLithotecQ0 quartz;
  bar.setVerbosity(verbosity);
  bar.setOptMaterial(quartz);
  bar.addSurface(a6);
  bar.addSurface(a2);
  bar.addSurface(a3);
  bar.addSurface(a4);
  bar.addSurface(a5);
  bar.addSurface(a1); // colored objects last to prevent black lines on top
  bar.setName("bar");


  // define the photon detector box
  //               7
  //             q3
  //            /|
  //           / |
  //          /  |
  //         /8  |
  //       q4    | 
  //       |     |
  //       |     |
  //       |     |
  //       |     |
  //       |     |
  //       | 5   |6
  //       q1----q2
  //
  // 
  XYZPoint q1(-100.0,  -5.0,-1300.0);
  XYZPoint q2(-100.0,  -5.0, -800.0);
  XYZPoint q3(-100.0, 595.0, -800.0);
  XYZPoint q4(-100.0, 595.0,-1300.0);

  XYZPoint q5( 100.0,  -5.0,-1300.0);
  XYZPoint q6( 100.0,  -5.0, -800.0);
  XYZPoint q7( 100.0, 595.0, -800.0);
  XYZPoint q8( 100.0, 595.0,-1300.0);


  // Declare flat surfaces with arbitrary number of points.
  PndDrcSurfPolyFlat b1,b2,b3,b4,b5,b6;

  b1.setVerbosity(verbosity);
  b1.addPoint(q1);
  b1.addPoint(q2);
  b1.addPoint(q3);
  b1.addPoint(q4);
  b1.setReflectivity(refl_none);
  b1.setName("bright");

  b2.setVerbosity(verbosity);
  b2.addPoint(q5);
  b2.addPoint(q6);
  b2.addPoint(q7);
  b2.addPoint(q8);
  b2.setReflectivity(refl_none);
  b2.setName("bleft");

  b3.setVerbosity(verbosity);
  b3.addPoint(q1);
  b3.addPoint(q5);
  b3.addPoint(q8);
  b3.addPoint(q4);
  b3.setReflectivity(refl_none);
  b3.setName("bup");

  b4.setVerbosity(verbosity);
  b4.addPoint(q2);
  b4.addPoint(q6);
  b4.addPoint(q7);
  b4.addPoint(q3);
  b4.setReflectivity(refl_none);
  b4.setName("bdown");

  b5.setVerbosity(verbosity);
  b5.addPoint(q4);
  b5.addPoint(q8);
  b5.addPoint(q7);
  b5.addPoint(q3);
  b5.setReflectivity(refl_none);
  b5.setName("btop");

  b6.setVerbosity(verbosity);
  b6.addPoint(q1);
  b6.addPoint(q2);
  b6.addPoint(q6);
  b6.addPoint(q5);
  b6.setReflectivity(refl_none);
  b6.setName("bbottom");

  // focussing mirror
  double    radius;
  double    half_width;
  radius     = (imirror==1) ? 400 : 500;
  half_width = (imirror==1) ? 210 : 300;
  XYZPoint  rim(0,-half_width,radius);
  XYZPoint  ref(rim);

  
  Transform3D trans1;
  trans1 = Transform3D(RotationY(pi));

  if (imirror==1)  trans1 *= Transform3D(RotationX(pi/180.0*(8)));
  if (imirror==2)  trans1 *= Transform3D(RotationX(pi/180.0*(0)));
  rim = trans1 * rim;
  XYZPoint pos1(0,-100,-1299);// here the position should be
  XYZVector shift = pos1-rim;
  trans1 = Transform3D(shift)*trans1; // multiply from the left!
  ref = trans1 * ref;

  PndDrcSurfPolySphere bmirror1a;
  bmirror1a.addPoint(XYZPoint(-half_width,-half_width,0));
  bmirror1a.addPoint(XYZPoint(-half_width,+half_width,0));
  bmirror1a.addPoint(XYZPoint(+half_width,+half_width,0));
  bmirror1a.addPoint(XYZPoint(+half_width,-half_width,0));
  bmirror1a.setRadius(radius);
  bmirror1a.addTransform(trans1);
  

  bmirror1a.setName("bmirror1a");
  bmirror1a.setPrintColor(4);
  bmirror1a.setReflectivity(refl_perfect);

  PndDrcSurfPolyPara bmirror1b;
  bmirror1b.addPoint(XYZPoint(-half_width,-half_width,0));
  bmirror1b.addPoint(XYZPoint(-half_width,+half_width,0));
  bmirror1b.addPoint(XYZPoint(+half_width,+half_width,0));
  bmirror1b.addPoint(XYZPoint(+half_width,-half_width,0));
  bmirror1b.setRadius(radius);
  bmirror1b.addTransform(trans1);
  bmirror1b.setName("bmirror1b");
  bmirror1b.setPrintColor(4);
  bmirror1b.setReflectivity(refl_perfect);


  // the second flat mirror 
  XYZPoint qq2(-100.0, -5.0,  -970);
  XYZPoint qq6( 100.0, -5.0,  -970);
  XYZPoint qq3(-100.0,150.0, -1000);
  XYZPoint qq7( 100.0,150.0, -1000);
  shift = XYZVector(0,0,60);
  qq2 += shift;
  qq6 += shift;
  qq3 += shift;
  qq7 += shift;

  PndDrcSurfPolyFlat bmirror2;  
  bmirror2.setVerbosity(verbosity);
  bmirror2.addPoint(qq2);
  bmirror2.addPoint(qq6);
  bmirror2.addPoint(qq7);
  bmirror2.addPoint(qq3);
  bmirror2.setName("bmirror2");
  bmirror2.setPrintColor(4);               
  bmirror2.setReflectivity(refl_perfect);

  // focussing mirror
  radius     = (imirror==1) ? 250 : 180;
  half_width = (imirror==1) ? 150 : 90;
  rim.SetXYZ(0,-half_width,radius);

  Transform3D trans2;
  if (imirror==1) trans2 = Transform3D(RotationX(pi/180.0*(15)));
  if (imirror==2) trans2 = Transform3D(RotationX(pi/180.0*(10)));
  rim = trans2 * rim;
  if (imirror==1) pos1.SetXYZ(0,140,-880);// here the position should be
  if (imirror==2) pos1.SetXYZ(0,320,-830);// here the position should be
  shift = pos1-rim;
  trans2 = Transform3D(shift) * trans2;

  PndDrcSurfPolySphere bmirror3a;
  bmirror3a.addPoint(XYZPoint(-half_width,-half_width,0));
  bmirror3a.addPoint(XYZPoint(-half_width,+half_width,0));
  bmirror3a.addPoint(XYZPoint(+half_width,+half_width,0));
  bmirror3a.addPoint(XYZPoint(+half_width,-half_width,0));
  bmirror3a.setRadius(radius);
  bmirror3a.addTransform(trans2);
  bmirror3a.setName("bmirror3a");
  bmirror3a.setPrintColor(4);
  bmirror3a.setReflectivity(refl_perfect);

  PndDrcSurfPolyPara bmirror3b;
  bmirror3b.addPoint(XYZPoint(-half_width,-half_width,0));
  bmirror3b.addPoint(XYZPoint(-half_width,+half_width,0));
  bmirror3b.addPoint(XYZPoint(+half_width,+half_width,0));
  bmirror3b.addPoint(XYZPoint(+half_width,-half_width,0));
  bmirror3b.setRadius(radius);
  bmirror3b.addTransform(trans2);
  bmirror3b.setName("bmirror3b");
  bmirror3b.setPrintColor(4);
  bmirror3b.setReflectivity(refl_perfect);


  
  double angle_screen = (imirror==1) ? 20 : 75;
  RotationX rot_screen(pi/180.0*angle_screen);
  XYZVector shift_screen;
  (imirror==1) ? shift_screen.SetXYZ(0,370,-950): shift_screen.SetXYZ(0,462,-910);
  half_width = (imirror==1) ? 100 : 100;
  PndDrcSurfPolyFlat bscreen;

  Transform3D trans3(rot_screen);
  trans3 = Transform3D(shift_screen)*trans3;

  bscreen.addPoint(XYZPoint(-half_width,-half_width,-10));
  bscreen.addPoint(XYZPoint(-half_width,+half_width,-10));
  bscreen.addPoint(XYZPoint(+half_width,+half_width,-10));
  bscreen.addPoint(XYZPoint(+half_width,-half_width,-10));

  //bscreen.addPoint(p5);
  //bscreen.addPoint(p6);
  //bscreen.addPoint(p7);
  //bscreen.addPoint(p8);



  bscreen.setName("bscreen");
  bscreen.setPrintColor(2);
  bscreen.setPixel();
  bscreen.addTransform(trans3);

  PndDrcOptVol box;
  box.setVerbosity(verbosity);
  box.setOptMaterial(quartz);
  box.addSurface(b1);
  box.addSurface(b2);
  box.addSurface(b3);
  box.addSurface(b4);
  box.addSurface(b5);
  box.addSurface(b6);
  if (imirror==1) box.addSurface(bmirror1a);
  if (imirror==2) box.addSurface(bmirror1b);
  if (imirror==1) box.addSurface(bmirror3a);
  if (imirror==2) box.addSurface(bmirror3b);
  box.addSurface(bscreen);
  box.setName("box");




  // Build a optical system consisting out of several volumes,
  // mirrors and screens.
  // This layer has the advantage, that a device consisting out 
  // of many equal subsystems
  // like a bar box, easily can be reproduced. 
  // See one of the forth comming test examples.

  PndDrcOptDevSys opt_system;
  opt_system.setName("total");
  opt_system.setVerbosity(verbosity);
  //opt_system.addDevice(cube);
  opt_system.addDevice(bar);
  opt_system.addDevice(box);

  // couple surface 1 of device 1 with surface 2 of device 2
  //                       dev1   dev2    surf1     surf2
  opt_system.coupleDevice("bar","box","aup",  "bdown");
  
  // The manager must be created as pointer. It is created as singleton, that is only 
  // one manager can exist per application.
  PndDrcOptDevManager* manager = new PndDrcOptDevManager();
  manager->setVerbosity(verbosity);
  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(-500,300,-1300,500,620,-500);"<<endl;
  //geo<<"    view->SetRange(-500,-300,-1300,500,700,1300);"<<endl;
  geo<<"    Int_t i;"<<endl;
  geo<<"    view->SetView(180,90,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;
  //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);
  //
  // the intention is to play around with routines.

  // create a list of photons in bar
  //exit(1);

  list<PndDrcPhoton> list_photon; // get list

  PndDrcPhoton ph;


  bool photons_exist = false;

  if (ioption==1)  
    {
      for (double xx=-9; xx<=9; xx+=1.0)
	{
	  for (double yy=-4; yy<=4; yy+=1.0)
	    {
	      ph.setPosition(XYZPoint(xx,yy,0));
	      ph.setDirection(XYZVector(0,0,-1));
	      ph.setWavelength(600);
	      list_photon.push_back(ph);
	    }
	}
      photons_exist = true;
      manager->setPhotonList(list_photon,"bar","total",0,0);
    }
  if (ioption==2)  
    {
      for (double xx=-9; xx<=9; xx+=3)
	{
	  for (double yy=-4; yy<=4; yy+=2)
	    {
	      //double xx=0;
	      //double yy=0;
	      ph.setPosition(XYZPoint(xx,yy,0));
	      ph.setDirection(XYZVector(0,0,-1));
	      ph.setWavelength(650);
	      list_photon.push_back(ph);
	      
	      double angle = pi/180.0 * 20.0;
	      ph.setDirection(XYZVector(0,sin(+angle),-cos(+angle)));
	      ph.setWavelength(521);
	      list_photon.push_back(ph);
	      ph.setDirection(XYZVector(0,sin(-angle),-cos(-angle)));
	      ph.setWavelength(521);
	      list_photon.push_back(ph);
	     
	      angle = pi/180.0 * 40.0;
	      ph.setDirection(XYZVector(0,sin(+angle),-cos(+angle)));
	      ph.setWavelength(451);
	      list_photon.push_back(ph);
	      ph.setDirection(XYZVector(0,sin(-angle),-cos(-angle)));
	      ph.setWavelength(451);
	      list_photon.push_back(ph);
	      
	    }
	}
      photons_exist = true;
      manager->setPhotonList(list_photon,"bar");
    }
  if (ioption==3)
    {
      XYZPoint pos(0,-10,1000);
      XYZVector dir(0,1,4); 
      double   beta = 0.69;
      double range = 500;
      photons_exist = manager->cerenkov(pos,dir,beta,100,range,400,405); // generate photons
    }


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

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


  // propagate writes to geo, that has finished, therefore, close geo
  geo<<"}"<<endl;
  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]\",200,-100,100);"<<endl;
  scr<<"    hgr->SetMarkerStyle(20);"<<endl;
  scr<<"    hgr->SetMarkerSize(0.2);"<<endl;
  scr<<"    hgr->SetMinimum(-100);"<<endl;
  scr<<"    hgr->SetMaximum(+100);"<<endl;
  scr<<"    hgr->Draw(\"POL\");"<<endl;
  int icnt_measured = 0;
  int icnt_flying   = 0;
  int icnt_lost     = 0;
  int icnt_absorbed = 0;

  Transform3D trans3_inv = trans3.Inverse();
  list<PndDrcPhoton>::iterator iph;
  for(iph=list_photon.begin(); iph != list_photon.end(); ++iph) 
    {
      if ((*iph).fate()==Drc::PhotMeasured)
	{
	  icnt_measured++;
	  XYZPoint pos = (*iph).position();
	  pos = trans3_inv * pos;

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


	  scr<<"    TMarker* t = new TMarker("<<x<<","<<y<<",20);"<<endl;
	  scr<<"    t->SetMarkerColor("
		<<(*iph).colorNumber((*iph).wavelength())
		<<");"<<endl;
	  scr<<"    t->SetMarkerSize(0.7);"<<endl;
	  scr<<"    t->Draw();"<<endl;
	}
      else if ((*iph).fate()==Drc::PhotFlying)   icnt_flying++; // should never happen.
      else if ((*iph).fate()==Drc::PhotAbsorbed) icnt_absorbed++;
      else                                       icnt_lost++;
    }
  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;

}