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

#include <map>
using std::map;

#include <utility>
using std::pair;

#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 "DrcPhoton.h"
#include "DrcOptReflAbs.h"
#include "DrcSurfAbs.h"
#include "DrcSurfPolyFlat.h"
#include "DrcOptReflSilver.h"
#include "DrcOptMatAbs.h"
#include "DrcOptMatLithotecQ0.h"
#include "DrcOptDev.h"
#include "DrcOptDevSys.h"
#include "DrcOptVol.h"
#include "DrcOptPixel.h"
#include "DrcOptMirror.h"
#include "DrcOptDevManager.h"


int main()
{

  // Example for a simple bar with screen (photon detection) and mirror.

  //int number_photons = 50;

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


  // create one segment of a disk and reprocuce the other.
  //
  //                                        y
  //                                        |
  //                                        |
  //                                        |
  //              8                  7      |
  //           p4-----------------p3        |
  //            \                /          0-----------x
  //             \              /          /
  //              \            /          / 
  //               \          /          /  
  //                \  5     / 6        z   
  //                 p1-----p2              
  //                                        
  
  const int nsegs      = 12;           // Number of segments
  double    dist_z1    = 200;//mm      // distance from beam axis (=z)
  double    dist_z2    = 700;//mm      // outer distance from beam.
  double    angle      = 2*3.1415926/nsegs; 


  XYZPoint p1(-dist_z1*tan(angle/2), dist_z1, 200); 
  XYZPoint p2(+dist_z1*tan(angle/2), dist_z1, 200); 
  XYZPoint p3(+dist_z2*tan(angle/2), dist_z2, 200); 
  XYZPoint p4(-dist_z2*tan(angle/2), dist_z2, 200); 
  XYZPoint p5(-dist_z1*tan(angle/2), dist_z1, 210); 
  XYZPoint p6(+dist_z1*tan(angle/2), dist_z1, 210); 
  XYZPoint p7(+dist_z2*tan(angle/2), dist_z2, 210); 
  XYZPoint p8(-dist_z2*tan(angle/2), dist_z2, 210); 

  DrcSurfPolyFlat pfront;
  //pfront.setVerbosity(4);
  pfront.addPoint(p1);  
  pfront.addPoint(p2);
  pfront.addPoint(p3);
  pfront.addPoint(p4);
  pfront.setName("pfront");

  DrcSurfPolyFlat pback(pfront); // construct from pfront
  pback.setVerbosity(verbosity);
  Transform3D trans1(XYZVector(0,0,+10));
  pback.addTransform(trans1);// and shift it
  pback.setName("pback");
  //DrcSurfPolyFlat pback;
  //pfront.setVerbosity(4);
  //pback.addPoint(p5);  
  //pback.addPoint(p6);
  //pback.addPoint(p7);
  //pback.addPoint(p8);
  //pback.setName("pback");

  DrcSurfPolyFlat pbottom;
  pbottom.addPoint(p1);  
  pbottom.addPoint(p2);
  pbottom.addPoint(p6);
  pbottom.addPoint(p5);
  pbottom.setName("pbottom");

  DrcSurfPolyFlat ptop;
  ptop.addPoint(p4);  
  ptop.addPoint(p3);
  ptop.addPoint(p7);
  ptop.addPoint(p8);
  ptop.setName("ptop");

  DrcSurfPolyFlat pleft;
  pleft.setVerbosity(verbosity);
  pleft.addPoint(p1);  
  pleft.addPoint(p5);
  pleft.addPoint(p8);
  pleft.addPoint(p4);
  pleft.setName("pleft");

  DrcSurfPolyFlat pright;
  pright.setVerbosity(verbosity);
  pright.addPoint(p2);  
  pright.addPoint(p6);
  pright.addPoint(p7);
  pright.addPoint(p3);
  pright.setName("pright");

  DrcOptVol psegment;
  DrcOptMatLithotecQ0 quartz;
  psegment.setOptMaterial(quartz);
  psegment.addSurface(pfront);
  psegment.addSurface(pback);
  psegment.addSurface(pbottom);
  psegment.addSurface(ptop);
  psegment.addSurface(pleft);
  psegment.addSurface(pright);
  psegment.setName("segment");
  psegment.setVerbosity(verbosity);

  DrcOptPixel pscreen;
  DrcSurfPolyFlat ps(ptop);
  ps.setName("screen_surf");
  ps.setPrintColor(2);         // red
  pscreen.addSurface(ps);
  pscreen.setName("screen");

  DrcSurfPolyFlat amf(pbottom);
  amf.setName("mirror_front");
  amf.setPrintColor(4);        // blue
  DrcSurfPolyFlat amb(pbottom);
  amb.setName("mirror_back");
  amb.setPrintColor(4);        // blue
  DrcOptMirror pmirror;
  pmirror.setFrontSurface(amf);
  pmirror.setBackSurface(amb);
  pmirror.setName("mirror");


  DrcOptDevSys psystem;
  psystem.addDevice(psegment);
  psystem.addDevice(pscreen);
  psystem.addDevice(pmirror);

  psystem.coupleDevice("segment","screen","ptop","screen_surf");
  psystem.coupleDevice("segment","mirror","pbottom","mirror_front");
  psystem.setName("sys_segment");

  DrcOptDevManager manager;
  // reproduce the system nsegs times (system = segment + mirror)
  DrcOptDevSys     psys[nsegs];



  for (int i=0; i<nsegs; i++)
  {
    psys[i] = psystem; // copy
    RotationZ   rotz(-angle*i); // clockwise otherwise pleft=pright 
    //                             coupl. below fails...
    cout<<" rotate "<<i<<" by "<<-angle*i*180/3.1415926<<endl;
    Transform3D trans2(rotz);
    //if (i==0) psys[i].setPrintColor(2);
    //if (i==1) psys[i].setPrintColor(3);
    psys[i].addTransform(trans2);
    manager.addDeviceSystem(psys[i],i);
  }


  // couple systems side at side
    for (int i=1; i<nsegs; i++)
      {
	manager.coupleDevice("segment","segment","pleft","pright",i,(i-1)); 
      }
    // connect last with first.
    manager.coupleDevice("segment","segment","pleft","pright",0,nsegs-1);

	
  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,-500,-500,500,500,500);"<<endl;
  geo<<"    Int_t i;"<<endl;
  geo<<"    view->SetView(0,90,90,i);"<<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.

  // photons

  // get volume and positions for intersection with particle.
  // A intersection with the volume produces on complex shapes 2n hits
  // with a volume name and a copy number.

  /*
  list<DrcPhoton> list_photon;
  DrcPhoton ph;
  ph.setPosition(XYZPoint(0,300,205));
  ph.setDirection(XYZVector(-2,1,0));
  ph.setWavelength(500);
  list_photon.push_back(ph);
  bool photons_exist = true;
  manager.setPhotonList(list_photon,"segment",0);
  */

  
    XYZPoint  pos(0,0,0);
    XYZVector dir(0,2,2);
    bool photons_exist = manager.cerenkov(pos,dir.Unit(),0.99); // generate photons
  
  //bool photons_exist = (photons_exist1 || photons_exist2 || 
  //		photons_exist3 || photons_exist4); 

  if (photons_exist)
    {
      fstream scr;
      scr.open("Screen.C",std::ios::out);
      scr<<"{"<<endl;
      scr<<"    TCanvas *c1 = new TCanvas(\"c1\"); "<<endl;
      scr<<"    TH1F *hgr = new TH1F(\"hgr\",\"a simple graph\",100,-100,100);"<<endl;
      scr<<"    hgr->SetMarkerStyle(20);"<<endl;
      scr<<"    hgr->SetMinimum(-100);"<<endl;
      scr<<"    hgr->SetMaximum(+100);"<<endl;
      scr<<"    hgr->Draw(\"POL\");"<<endl;



      int icnt_measured = 0;
      int icnt_absorbed = 0;
      int icnt_lost     = 0;

      manager.propagate();                           // propagate photons

      list<DrcPhoton> list_photon = manager.photonList(); // get list

      list<DrcPhoton>::iterator iph;
      for (iph = list_photon.begin(); iph != list_photon.end(); ++iph)
	{
	  //cout<<(*iph).fate()<<endl;
	  if ((*iph).fate()==Drc::PhotMeasured)
	    {
	      icnt_measured++;
	      double x     = (*iph).direction().X();
	      double y     = (*iph).direction().Y();
	      double z     = (*iph).direction().Z();
	      double theta = atan2(sqrt(x*x+y*y),fabs(z));
	      double phi   = atan2(y,x);
	      double r     = theta;
	      double xx = r*cos(phi)*180/3.1415;
	      double yy = r*sin(phi)*180/3.1415;
	      scr<<"    TMarker* t = new TMarker("<<xx<<","<<yy<<",20);"<<endl;
	      scr<<"    t->SetMarkerColor("
		 <<(*iph).colorNumber((*iph).wavelength())
		 <<");"<<endl;
	      scr<<"    t->SetMarkerSize(0.7);"<<endl;
	      scr<<"    t->Draw();"<<endl;
	    }
	  if ((*iph).fate()==Drc::PhotAbsorbed) icnt_absorbed++;
	  if ((*iph).fate()==Drc::PhotLost)     icnt_lost++;

	  
	}
      scr<<"}"<<endl;
      scr.close();
      int icnt = icnt_measured+icnt_absorbed+icnt_lost;
      cout<<" Total photons: "<<icnt<<endl;
      cout<<"      measured: "<<icnt_measured<<endl;
      cout<<"      absorbed: "<<icnt_absorbed<<endl;
      cout<<"          lost: "<<icnt_lost<<endl;
    }



  geo<<"}"<<endl;
  geo.close();
  
  
  return 0;
}