#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 "TRint.h"
#include "TFile.h"
#include "TVector3.h"
#include "TRandom.h"
#include "TRotation.h"
#include "TNtuple.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 "PndDrcOptReflSilver.h"
#include "PndDrcOptReflGray.h"
#include "PndDrcOptMatAbs.h"
#include "PndDrcOptMatLithotecQ0.h"
#include "PndDrcOptDev.h"
#include "PndDrcOptDevSys.h"
#include "PndDrcOptVol.h"
#include "PndDrcOptDevManager.h"


int main(int argc, char *argv[])
{

  //TRint *gROOT = new TRint("ROOT example", &argc, argv);

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

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


  // 8 points define a bar

  XYZPoint p1(-200.0, +5.0, 0);    //         p5----------p8
  XYZPoint p2(-200.0, -5.0, 0);    //        /|          /|
  XYZPoint p3(+200.0, -5.0, 0);    //       / |         / |
  XYZPoint p4(+200.0, +5.0, 0);    //      /  p6-------/--p7
  //                                     /  /        /  /
  //                                    /  /        /  /
  //                                   /  /        /  /
  //                                  /  /        /  /
  XYZPoint p5(-200.0, +5.0, -1000);  // p1---------p4  /
  XYZPoint p6(-200.0, -5.0, -1000);  // | /         | /
  XYZPoint p7(+200.0, -5.0, -1000);  // |/          |/
  XYZPoint p8(+200.0, +5.0, -1000);  // 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.

  PndDrcOptReflSilver refl;

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

  a1.setVerbosity(verbosity);
  a1.setReflectivity(refl);  // mirror
  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.setPixel();         // screen
  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(a1);
  bar.addSurface(a2);
  bar.addSurface(a3);
  bar.addSurface(a4);
  bar.addSurface(a5);
  bar.addSurface(a6);
  bar.setName("bar");


  /*
  PndDrcSurfPolyFlat bmf,bmb;
  PndDrcOptReflGray refl_semi;
  refl_semi.setReflProb(0.7);
  bmf = a6;
  bmf.setReflectivity(refl_semi);
  bmf.setName("bmirror_front");
  bmf.setPrintColor(4);
  bmb = a6;
  bmb.setName("bmirror_back");
  bmb.setPixel();               // screen
  bmb.setPrintColor(2);

  PndDrcOptVol mirror1;
  mirror1.setOptMaterial(quartz);
  mirror1.addSurface(bmf);
  mirror1.addSurface(bmb);
  mirror1.setName("mirror1");
  mirror1.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.

  PndDrcOptDevSys opt_system;
  opt_system.setVerbosity(verbosity);
  opt_system.addDevice(bar);
  //opt_system.addDevice(mirror1);



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

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

  bool l_vis=false;
  
  if (l_vis)
    {
      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(-900,-900,-900,900,900,900);"<<endl;
      geo<<"    Int_t i;"<<endl;
      geo<<"    view->SetView(0,90,90,i);"<<endl;
      manager->print(geo);
    }
  
  // create a list of photons in bar

  XYZPoint pos(0,-10,-50);
  XYZVector dir(0,0.3,-1); 


  TFile *f = new TFile("test_top.root","RECREATE");
  TNtuple *ntuple = new TNtuple("ntuple","test_top data","x:y:z:t:m");  

  bool photons_exist=false;
  
  bool l_dispersion=false;

  int nph = 5000;

  
  // pions: 
  double   beta = 1.0/sqrt(1.0+0.140*0.140);
  if (l_dispersion)
    {
      photons_exist = manager->cerenkov(pos,dir,beta,nph,1000,400,600); // gen photons
    }
  else
    {
      photons_exist = manager->cerenkov(pos,dir,beta,nph,1000,400,400); // gen photons
    }
  if (photons_exist) manager->propagate();                   // propagate photons
  list<PndDrcPhoton> list_photon = manager->photonList(); // get list
  list<PndDrcPhoton>::iterator iph;
  for(iph=list_photon.begin(); iph != list_photon.end(); ++iph) 
    {
      if ((*iph).fate()==Drc::PhotMeasured)
	{
	  double x     = (*iph).position().X();
	  double y     = (*iph).position().Y();
	  double z     = (*iph).position().Z();
	  double t     = (*iph).time();
	  double m=140.0;
	  ntuple->Fill(x,y,z,t,m);
	}
    }
  manager->clearPhotonList();
  // kaons
  beta = 1.0/sqrt(1.0+0.500*0.500);
  
  
  if (l_dispersion)
    {
      photons_exist = manager->cerenkov(pos,dir,beta,nph,1000,400,600); // gen photons
    }
  else
    {
      photons_exist = manager->cerenkov(pos,dir,beta,nph,1000,400,400); // gen photons
    }
  
  if (photons_exist) manager->propagate();                   // propagate photons
  list_photon.clear();
  list_photon = manager->photonList(); // get list
  for(iph=list_photon.begin(); iph != list_photon.end(); ++iph) 
    {
      if ((*iph).fate()==Drc::PhotMeasured)
	{
	  double x     = (*iph).position().X();
	  double y     = (*iph).position().Y();
	  double z     = (*iph).position().Z();
	  double t     = (*iph).time();
	  double m=500.0;
	  ntuple->Fill(x,y,z,t,m);
	}
    }



  f->Write();

  if (l_vis)
    {    
      // propagate writes to geo, that has finished, therefore, close geo
      geo<<"}"<<endl;
      geo.close();
    }
  
  delete manager;


 


  return EXIT_SUCCESS;

}