#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 "PndDrcOptReflSurf.h"
#include "PndDrcSurfAbs.h"
#include "PndDrcSurfPolyFlat.h"
#include "PndDrcOptReflSilver.h"
#include "PndDrcOptMatAbs.h"
#include "PndDrcOptMatLithotecQ0.h"
#include "PndDrcOptDev.h"
#include "PndDrcOptDevSys.h"
#include "PndDrcOptVol.h"
#include "PndDrcOptDevManager.h"


int main()
{

  PndDrcOptReflSurf refl_surf;
  
  string filename;
  double prob;
  
  cout<< "Give filename and probability" <<endl;
  cin>>filename;
  cin>>prob;
  


  refl_surf.setProb(prob);
  
  cout<<" probability: "<<refl_surf.prob()<<endl;
  


  TRandom ran;


  // 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(-17.0, +8.0, +400);    //         p5----------p8
  XYZPoint p2(-17.0, -8.0, +400);    //        /|          /|
  XYZPoint p3(+17.0, -8.0, +400);    //       / |         / |
  XYZPoint p4(+17.0, +8.0, +400);    //      /  p6-------/--p7
  //                                     /  /        /  /
  //                                    /  /        /  /
  //                                   /  /        /  /
  //                                  /  /        /  /
  XYZPoint p5(-17.0, +8.0, -400);  // p1---------p4  /
  XYZPoint p6(-17.0, -8.0, -400);  // | /         | /
  XYZPoint p7(+17.0, -8.0, -400);  // |/          |/
  XYZPoint p8(+17.0, +8.0, -400);  // 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");

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


  a6.setPixel();
  a1.setPixel();

  a2.setReflectivity(refl_surf);
  a3.setReflectivity(refl_surf);
  a4.setReflectivity(refl_surf);
  a5.setReflectivity(refl_surf);
  
  // 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");

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


  
  // 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(-100,-100,-100,100,100,100);"<<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




  fstream out;
  out.open(filename.c_str(),std::ios::out);

  for (double zz=-390; zz<=390; zz+=10)
    {

      XYZPoint pos(0,-9,zz);
      double nevents=1000;
      double pi=3.1415926535;



      bool photons_exist = false;
      int icnt_coinc    = 0;
      int icnt_measured = 0;
      int icnt_flying   = 0;
      int icnt_lost     = 0;
      int icnt_absorbed = 0;

      for (int ievent=0; ievent<nevents; ievent++)
	{

	  double costh = ran.Uniform(0,1);// 90-0 deg
	  double phi   = ran.Uniform(0,2*pi);
	  double sinth = sqrt(1-costh*costh);
	  // emission in half sphere
      
	  XYZVector dir(sinth*sin(phi),costh,sinth*cos(phi)); 
	  double    beta = 0.8;
      
	  if (manager->cerenkov(pos,dir,beta,300,200,400,600)) photons_exist = true; 
  
	  if (photons_exist)
	    {
	      manager->propagate();                           // propagate photons
	    }

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



	  bool det1 = false;
	  bool det2 = false;
	  double time1=0;
	  double time2=0;
	  
	  list<PndDrcPhoton>::iterator iph;
	  for(iph=list_photon.begin(); iph != list_photon.end(); ++iph) 
	    {
	      //if ((*iph).fate()==PndDrc::PhotMeasured) cout<<" Measured"<<endl;
	      //if ((*iph).fate()==PndDrc::PhotFlying)   cout<<" Flying"<<endl;
	      //if ((*iph).fate()==PndDrc::PhotAbsorbed) cout<<" Absorbed"<<endl;
	      //if ((*iph).fate()==PndDrc::PhotLost)     cout<<" Lost"<<endl;
	      


	      if ((*iph).fate()==Drc::PhotMeasured)
		{
		  icnt_measured++;
		  //double x     = (*iph).position().X();
		  //double y     = (*iph).position().Y();
		  double z     = (*iph).position().Z();
		  double t     = (*iph).time();
		  
		  if (fabs(-400-z)<0.1) {det1=true;time1=t;}
		  
		  if (fabs(+400-z)<0.1) {det2=true;time2=t;}

		  //cout<<z<<" "<<(*iph).time()<<endl;
		}
	      else if ((*iph).fate()==Drc::PhotFlying)   icnt_flying++; // should never happen.
	      else if ((*iph).fate()==Drc::PhotAbsorbed) icnt_absorbed++;
	      else                                       icnt_lost++;
	    }

	  if (det1 && det2 /*&& (time2 > time1)*/)
	    {
	      icnt_coinc++;
	      //cout<<" coincident hit of both detectors"<<endl;
	    }
	  else
	    {
	      //cout<<" no coinc."<<endl;
	    }
	  // cout<<" event = "<<ievent<<endl;
	}
      out<<zz<<" "<<icnt_coinc<<endl;
      cout<<zz<<" "<<icnt_coinc<<endl;
    }
  //int icnt = icnt_measured+icnt_flying+icnt_lost+icnt_absorbed;
  //cout<<" generated photons: "<<icnt<<endl;
  //cout<<" measured  photons: "<<icnt_measured<<endl;
  //cout<<" coinc.    photons: "<<icnt_coinc<<endl;
  //cout<<" absorbed  photons: "<<icnt_absorbed<<endl;
  //cout<<" lost      photons: "<<icnt_lost<<endl;


  // propagate writes to geo, that has finished, therefore, close geo
  geo<<"}"<<endl;
  geo.close();

  delete manager;

  return EXIT_SUCCESS;

}