#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 "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 verbosity = 3; // 0=quiet, 1=constructors,2=member,3=functionality
  //                    4=photons, 5=everything.


  // 8 points define a bar

  XYZPoint p1(-10.0, +5.0, 0);    //         p5----------p8
  XYZPoint p2(-10.0, -5.0, 0);    //        /|          /|
  XYZPoint p3(+10.0, -5.0, 0);    //       / |         / |
  XYZPoint p4(+10.0, +5.0, 0);    //      /  p6-------/--p7
  //                                     /  /        /  /
  //                                    /  /        /  /
  //                                   /  /        /  /
  //                                  /  /        /  /
  XYZPoint p5(-10.0, +5.0, -40);  // p1---------p4  /
  XYZPoint p6(-10.0, -5.0, -40);  // | /         | /
  XYZPoint p7(+10.0, -5.0, -40);  // |/          |/
  XYZPoint p8(+10.0, +5.0, -40);  // 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.
  DrcSurfPolyFlat a1,a2,a3,a4,a5,a6;

  a1.setVerbosity(verbosity);
  a1.addPoint(p1);
  a1.addPoint(p2);
  a1.addPoint(p3);
  a1.addPoint(p4);
  a1.setName("adown");

  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");

  DrcSurfPolyFlat as;
  as = a6;
  as.setName("ascreen");
  as.setPrintColor(2);                     // root red

  // A mirror has two sides, since behind a dichroic mirror a volume or screen
  // can follow.

  DrcSurfPolyFlat amf,amb;
  DrcOptReflSilver refl;
  amf = a1;
  amf.setReflectivity(refl);
  amf.setName("amirror_front");
  amf.setPrintColor(4);
  amb = a1;
  amb.setName("amirror_back");
  amb.setPrintColor(4);


  // create a volume consiting of surfaces
  // create a material the bar will consist of
  DrcOptVol bar;
  DrcOptMatLithotecQ0 quartz;
  bar.setVerbosity(verbosity);
  bar.setOptMaterial(quartz);
  bar.addSurface(a1);
  bar.addSurface(a2);
  bar.addSurface(a3);
  bar.addSurface(a4);
  bar.addSurface(a5);
  bar.addSurface(a6);
  bar.setName("bar");

  // create a screen (simplified photon detector) where photons end and 
  // get the status Drc::Measured (no further propagation). For this you need 
  // one single plane.

  DrcOptPixel screen;
  screen.setName("screen");
  screen.addSurface(as); 
  screen.setVerbosity(verbosity);

  DrcOptMirror mirror;
  mirror.setFrontSurface(amf);
  mirror.setBackSurface(amb);
  mirror.setName("mirror");
  mirror.setVerbosity(verbosity);


  // 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.

  DrcOptDevSys opt_system;
  opt_system.setVerbosity(verbosity);
  opt_system.addDevice(bar);
  opt_system.addDevice(screen);
  opt_system.addDevice(mirror);


  // couple surface 1 of device 1 with surface 2 of device 2
  //                       dev1   dev2    surf1     surf2
  opt_system.coupleDevice("bar","screen","aup",  "ascreen");
  opt_system.coupleDevice("bar","mirror","adown","amirror_front");

  
  // The manager must be created as pointer. It is created as singleton, that is only 
  // one manager can exist per application.
  DrcOptDevManager* manager = new DrcOptDevManager();
  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(-50,-50,-50,50,50,50);"<<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.


  // create a list of photons in bar

  XYZPoint pos(0,-10,-20);
  XYZVector dir(0,1,0); 
  double   beta = 0.99;

  bool photons_exist = manager->cerenkov(pos,dir,beta); // generate photons

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

  list<DrcPhoton> 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\",\"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_flying   = 0;
  int icnt_lost     = 0;
  int icnt_absorbed = 0;

  list<DrcPhoton>::iterator iph;
  for(iph=list_photon.begin(); iph != list_photon.end(); ++iph) 
    {
      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;
	}
      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;

}