#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 "DrcSurfQuadFlatDiff.h"
#include "DrcSurfPolySphere.h"
#include "DrcSurfPolyPara.h"
#include "DrcSurfPolyAsphere.h"
#include "DrcOptReflSilver.h"
#include "DrcOptReflNone.h"
#include "DrcOptReflPerfect.h"
#include "DrcOptMatAbs.h"
#include "DrcOptMatLithotecQ0.h"
#include "DrcOptMatNLAK33A.h"
#include "DrcOptMatLLF1.h"
#include "DrcOptMatVacuum.h"
#include "DrcOptDev.h"
#include "DrcOptDevSys.h"
#include "DrcOptVol.h"
#include "DrcOptPixel.h"
#include "DrcOptMirror.h"
#include "DrcOptDevManager.h"
#include "DrcUtil.h"


int main(int argc, char *argv[])
{
  double conical_const = 0;


  //cout<<" constant = ";
  //cin>>conical_const;

  //double conical_const = -1; // parabola
  //double conical_const = 1; // ellipse
  //double conical_const = -1; // hyperbola

  //int ioption = 1; // parallel beam (0 deg)
  int ioption = 2; // 0,  +- 20,  +-40 deg
  //int ioption = 3; // C-cone
  //int ioption = 4; // random in front of lens.
  //int ioption = 5; // single photon for debugging.
  //int ioption = 6; // grid 5 deg

  double pi=3.1415926535;

  DrcOptReflPerfect refl_perfect;
  DrcOptReflNone    refl_none;
  double lambda=500;

  double slab_width  = 17;
  double slab_height = 17;


  XYZPoint q0(-slab_width/2, +slab_height/2, 0.0);
  XYZPoint q1(-slab_width/2, -slab_height/2, 0.0);
  XYZPoint q2(+slab_width/2, -slab_height/2, 0.0);
  XYZPoint q3(+slab_width/2, +slab_height/2, 0.0);

  DrcOptMatLithotecQ0 quartz;
  DrcOptMatNLAK33A    nlak33a;
  DrcOptMatVacuum     vacuum;
  DrcOptMatLLF1       llf1;

  DrcSurfPolySphere lens_sphere;
  lens_sphere.setVerbosity(0);
  //lens_sphere.setConicalConstant(conical_const);
  lens_sphere.addPoint(q0);
  lens_sphere.addPoint(q1);
  lens_sphere.addPoint(q2);
  lens_sphere.addPoint(q3);
  // screen 101 mm from b0
  double f,f1,f2,radius1,radius2;
  double n2 = llf1.refIndex(lambda);
  double n1 = nlak33a.refIndex(lambda);

  cout<<" n1 = "<<n1<<endl;
  cout<<" n2 = "<<n2<<endl;

  f  = 200;
  //cout<<" Give focal length : ";
  //cin>>f;

  double viol = 0.1;             // Petzval violation
  f1 = f*(1.0+n2/(-n1));
  f2 = f*(1.0+n1/(-n2));

  cout<<" petzval = "<<f1*n1+f2*n2<<endl;
  cout<<" f  = "<<1.0/(1.0/f1+1.0/f2)<<endl;

  cout<<" f1 = "<<f1<<endl;
  cout<<" f2 = "<<f2<<endl;

  radius1 = -f1*(n2-n1)/n2;
  radius2 = -f2*(n1-n2)/n1;

  radius1 = f1*(n1-1);

  radius2 = 1.0/(f2*(n2-1)) + 1.0/radius1;
  radius2 = 1.0/radius2; 
  
  //radius1 = -19.5;
  //radius2 = 50;

  //cout<<" Give R1 (<0) : "<<endl;
  //radius1 = atof(argv[1]);
  //cout<<" Give R2 (>0) : "<<endl;
  //radius2 = atof(argv[2]);

  radius1 =-16.1;
  radius2 = 34.88;

  //radius1 = -999999;
  //radius2 = 1.0/(288*(n2-1));
  //radius2 = 1.0/radius2; 

  cout<<" r1 = "<<radius1<<endl;
  cout<<" r2 = "<<radius2<<endl;


  lens_sphere.setRadius(fabs(radius1)); // f = R/(n-1) = 100 mm
  //lens_sphere.setVerbosity(4);
  lens_sphere.setPrintColor(4);
  lens_sphere.setName("lens_sphere");

  RotationY rot(pi);
  Rotation3D rot_lens(rot); 
  Transform3D trans_rot_lens(rot_lens);
  RotationY rot2(0);
  Rotation3D rot_lens2(rot2); 
  Transform3D trans_rot_lens2(rot_lens2);
  lens_sphere.addTransform(trans_rot_lens);

  double dist = lens_sphere.centerPoint().Z();

  XYZVector shift_lens(0,0,-dist+0.2);

  
  Transform3D trans_shift_lens = Transform3D(shift_lens);
  lens_sphere.addTransform(trans_shift_lens);

  //lens_sphere.shift(XYZVector(0,0,-(dist-2)));






  // Declare flat surfaces with arbitrary number of points.
  q0.SetZ(0); // later be rotated by 180 deg.
  q1.SetZ(0);
  q2.SetZ(0);
  q3.SetZ(0);
  q0 = trans_rot_lens * q0;
  q1 = trans_rot_lens * q1;
  q2 = trans_rot_lens * q2;
  q3 = trans_rot_lens * q3;

  DrcSurfPolyFlat lens_base;
  lens_base.setVerbosity(0);
  lens_base.addPoint(q0);
  lens_base.addPoint(q1);
  lens_base.addPoint(q2);
  lens_base.addPoint(q3);
  lens_base.setPrintColor(4);
  lens_base.setName("lens_base");

  

  DrcSurfQuadFlatDiff lens_side1;
  lens_side1.addSurface(lens_base,  q0,q1);
  lens_side1.addSurface(lens_sphere,
			lens_sphere.limitingPoint(0),lens_sphere.limitingPoint(1));
  lens_side1.setPrintColor(4);
  lens_side1.setName("lens_side1");
  DrcSurfQuadFlatDiff lens_side2;
  lens_side2.setVerbosity(0);
  lens_side2.addSurface(lens_base,  q1,q2);
  lens_side2.addSurface(lens_sphere,
			lens_sphere.limitingPoint(1),lens_sphere.limitingPoint(2));
  lens_side2.setPrintColor(4);
  lens_side2.setName("lens_side2");
  DrcSurfQuadFlatDiff lens_side3;
  lens_side3.addSurface(lens_base,  q2,q3);
  lens_side3.addSurface(lens_sphere,
			lens_sphere.limitingPoint(2),lens_sphere.limitingPoint(3));
  lens_side3.setPrintColor(4);
  lens_side3.setName("lens_side3");
  DrcSurfQuadFlatDiff lens_side4;
  lens_side4.addSurface(lens_base,  q3,q0);
  lens_side4.addSurface(lens_sphere,
			lens_sphere.limitingPoint(3),lens_sphere.limitingPoint(0));
  lens_side4.setPrintColor(4);
  lens_side4.setName("lens_side4");


  lens_side1.setReflectivity(refl_none);
  lens_side2.setReflectivity(refl_none);
  lens_side3.setReflectivity(refl_none);
  lens_side4.setReflectivity(refl_none);
  

  DrcOptVol lens;
  lens.setVerbosity(0);
  lens.addSurface(lens_base);
  lens.addSurface(lens_sphere);
  lens.addSurface(lens_side1);
  lens.addSurface(lens_side2);
  lens.addSurface(lens_side3);
  lens.addSurface(lens_side4);
  lens.setOptMaterial(nlak33a);
  //lens.setOptMaterial(quartz);
  lens.setName("lens");

 

  // create a second lense
  DrcSurfPolySphere lens2_sphere1(lens_sphere);
  //lens2_sphere1.setConicalConstant(conical_const);
  lens2_sphere1.setVerbosity(0);
  lens2_sphere1.setPrintColor(4);
  lens2_sphere1.setName("lens2_sphere1");



  DrcSurfPolySphere lens2_sphere2;
  //lens2_sphere2.setConicalConstant(conical_const);
  lens2_sphere2.setVerbosity(0);
  lens2_sphere2.setPrintColor(3);

  q0 = trans_rot_lens2 * q0;
  q1 = trans_rot_lens2 * q1;
  q2 = trans_rot_lens2 * q2;
  q3 = trans_rot_lens2 * q3;
  q0.SetZ(0); // 1 mm appart from sphere edge
  q1.SetZ(0);
  q2.SetZ(0);
  q3.SetZ(0);

  lens2_sphere2.addPoint(q0);
  lens2_sphere2.addPoint(q1);
  lens2_sphere2.addPoint(q2);
  lens2_sphere2.addPoint(q3);

  lens2_sphere2.setRadius(fabs(radius2));// f = R/(n-1) = 100 mm

  lens2_sphere2.setPrintColor(3);
  lens2_sphere2.setName("lens2_sphere2");
  lens2_sphere2.addTransform(trans_rot_lens2);


  dist = lens_sphere.limitingPoint(0).Z();
  double dist2 = lens2_sphere2.limitingPoint(0).Z();
  XYZVector shift_lens2(0,0,dist-dist2+0.5);
  lens2_sphere2.addTransform(Transform3D(XYZVector(shift_lens2)));

  DrcSurfQuadFlatDiff lens2_side1;
  lens2_side1.addSurface(lens2_sphere1,
			lens2_sphere1.limitingPoint(0),lens2_sphere1.limitingPoint(1));
  lens2_side1.addSurface(lens2_sphere2,
			lens2_sphere2.limitingPoint(0),lens2_sphere2.limitingPoint(1));
  lens2_side1.setPrintColor(3);
  lens2_side1.setName("lens2_side1");
  DrcSurfQuadFlatDiff lens2_side2;
  lens2_side2.addSurface(lens2_sphere2,
			lens2_sphere1.limitingPoint(1),lens2_sphere1.limitingPoint(2));
  lens2_side2.addSurface(lens2_sphere2,
			lens2_sphere2.limitingPoint(1),lens2_sphere2.limitingPoint(2));
  lens2_side2.setPrintColor(3);
  lens2_side2.setName("lens2_side2");
  DrcSurfQuadFlatDiff lens2_side3;
  lens2_side3.addSurface(lens2_sphere1,
			lens2_sphere1.limitingPoint(2),lens2_sphere1.limitingPoint(3));
  lens2_side3.addSurface(lens2_sphere2,
			lens2_sphere2.limitingPoint(2),lens2_sphere2.limitingPoint(3));
  lens2_side3.setPrintColor(3);
  lens2_side3.setName("lens2_side3");
  DrcSurfQuadFlatDiff lens2_side4;
  lens2_side4.setVerbosity(0);
  lens2_side4.addSurface(lens2_sphere1,
			lens2_sphere1.limitingPoint(3),lens2_sphere1.limitingPoint(0));
  lens2_side4.addSurface(lens2_sphere2,
			lens2_sphere2.limitingPoint(3),lens2_sphere2.limitingPoint(0));
  lens2_side4.setPrintColor(3);
  lens2_side4.setName("lens2_side4");


  lens2_side1.setReflectivity(refl_none);
  lens2_side2.setReflectivity(refl_none);
  lens2_side3.setReflectivity(refl_none);
  lens2_side4.setReflectivity(refl_none);

  DrcOptVol lens2;
  lens2.setVerbosity(0);
  lens2.addSurface(lens2_sphere1);
  lens2.addSurface(lens2_sphere2);
  lens2.addSurface(lens2_side1);
  lens2.addSurface(lens2_side2);
  lens2.addSurface(lens2_side3);
  lens2.addSurface(lens2_side4);
  lens2.setOptMaterial(llf1);
  //lens2.setOptMaterial(quartz);
  //lens2.setOptMaterial(vacuum);
  lens2.setName("lens2");

 
  double thick1 = lens_sphere.centerPoint().Z();
  cout<<" thickness lens1 = "<<thick1<<endl;
  double thick2 = lens2_sphere2.centerPoint().Z()-lens2_sphere1.centerPoint().Z();
  cout<<" thickness lens2 = "<<thick2<<endl;
  // create a box with a hole for the lense

  // recalculate focal length       >0
  //double ff2 = 1.0/((n2-1)*(1.0/(-radius1)-1.0/(-radius2)
  //	    +(n1-1)*thick1/(n1*(-radius1)*(-radius2))));
//cout<<" ff2 = "<<ff2<<endl;
  

  // hole:
  DrcSurfQuadFlatDiff hole_side1(lens_side1);
  hole_side1.setPrintColor(1);
  hole_side1.setName("hole_side1");
  DrcSurfQuadFlatDiff hole_side2(lens_side2);
  hole_side2.setPrintColor(1);
  hole_side2.setName("hole_side2");
  DrcSurfQuadFlatDiff hole_side3(lens_side3);
  hole_side3.setPrintColor(1);
  hole_side3.setName("hole_side3");
  DrcSurfQuadFlatDiff hole_side4(lens_side4);
  hole_side4.setVerbosity(0);
  hole_side4.setPrintColor(1);
  hole_side4.setName("hole_side4");
  hole_side4.setNameSurface1("hole_base");
  hole_side4.setNameSurface2("hole_sphere");

  DrcSurfPolySphere hole_sphere(lens_sphere);
  hole_sphere.setVerbosity(0);
  hole_sphere.setPrintColor(1);
  hole_sphere.setName("hole_sphere");
  // hole2:
  DrcSurfQuadFlatDiff hole2_side1(lens2_side1);
  hole2_side1.setPrintColor(1);
  hole2_side1.setName("hole2_side1");
  DrcSurfQuadFlatDiff hole2_side2(lens2_side2);
  hole2_side2.setPrintColor(1);
  hole2_side2.setName("hole2_side2");
  DrcSurfQuadFlatDiff hole2_side3(lens2_side3);
  hole2_side3.setPrintColor(1);
  hole2_side3.setName("hole2_side3");
  DrcSurfQuadFlatDiff hole2_side4(lens2_side4);
  hole2_side4.setPrintColor(1);
  hole2_side4.setName("hole2_side4");
  DrcSurfPolySphere hole2_sphere2(lens2_sphere2);
  hole2_sphere2.setVerbosity(0);
  hole2_sphere2.setPrintColor(1);
  hole2_sphere2.setName("hole2_sphere2");


  XYZPoint p1(-300.0, +300.0, 0);  //            p5----------p8
  XYZPoint p2(-300.0, -300.0, 0);  //           /|          /|
  XYZPoint p3(+300.0, -300.0, 0);  //          / |         / |
  XYZPoint p4(+300.0, +300.0, 0);  //         /  p6-------/--p7
  //                                              /  /        /  /
  //                                             /  /        /  /
  //                                            /  /        /  /
  //                                           /  /        /  /
  XYZPoint p5(-300.0, +300.0, +288); // p1-s1--s2---p4 /
  XYZPoint p6(-300.0, -300.0, +288); // | /s5--s6   | /
  XYZPoint p7(+300.0, -300.0, +288); // |/ s7--s8   |/
  XYZPoint p8(+300.0, +300.0, +288); // p2-s3--s4--p3 
  //cout<<" f="<<f<<endl;
  //XYZPoint p5(-300.0, +300.0, f); // p1-s1--s2---p4 /
  //XYZPoint p6(-300.0, -300.0, f); // | /s5--s6   | /
  //XYZPoint p7(+300.0, -300.0, f); // |/ s7--s8   |/
  //XYZPoint p8(+300.0, +300.0, f); // p2-s3--s4--p3 
  XYZPoint s1(-slab_width/2, +300, 0);
  XYZPoint s2(+slab_width/2, +300, 0);
  XYZPoint s3(-slab_width/2, -300, 0);
  XYZPoint s4(+slab_width/2, -300, 0);
  XYZPoint s5(-slab_width/2, +slab_height/2, 0);
  XYZPoint s6(+slab_width/2, +slab_height/2, 0);
  XYZPoint s7(-slab_width/2, -slab_height/2, 0);
  XYZPoint s8(+slab_width/2, -slab_height/2, 0);

  DrcSurfPolyFlat box_side1;
  box_side1.addPoint(p1);
  box_side1.addPoint(p2);
  box_side1.addPoint(s3);
  box_side1.addPoint(s1);
  box_side1.setName("box_side1");
  DrcSurfPolyFlat box_side2;
  box_side2.addPoint(s2);
  box_side2.addPoint(p4);
  box_side2.addPoint(p3);
  box_side2.addPoint(s4);
  box_side2.setName("box_side2");
  DrcSurfPolyFlat box_side3;
  box_side3.addPoint(s1);
  box_side3.addPoint(s2);
  box_side3.addPoint(s6);
  box_side3.addPoint(s5);
  box_side3.setName("box_side3");
  DrcSurfPolyFlat box_side4;
  box_side4.addPoint(s7);
  box_side4.addPoint(s8);
  box_side4.addPoint(s4);
  box_side4.addPoint(s3);
  box_side4.setName("box_side4");

  DrcSurfPolyFlat box_side5;
  box_side5.addPoint(p1);
  box_side5.addPoint(p2);
  box_side5.addPoint(p6);
  box_side5.addPoint(p5);
  box_side5.setReflectivity(refl_none);
  box_side5.setName("box_side5");
  DrcSurfPolyFlat box_side6;
  box_side6.addPoint(p2);
  box_side6.addPoint(p3);
  box_side6.addPoint(p7);
  box_side6.addPoint(p6);
  box_side6.setReflectivity(refl_none);
  box_side6.setName("box_side6");
  DrcSurfPolyFlat box_side7;
  box_side7.addPoint(p3);
  box_side7.addPoint(p4);
  box_side7.addPoint(p8);
  box_side7.addPoint(p7);
  box_side7.setReflectivity(refl_none);
  box_side7.setName("box_side7");
  DrcSurfPolyFlat box_side8;
  box_side8.addPoint(p1);
  box_side8.addPoint(p4);
  box_side8.addPoint(p8);
  box_side8.addPoint(p5);
  box_side8.setReflectivity(refl_none);
  box_side8.setName("box_side8");
  DrcSurfPolyFlat box_side9;
  box_side9.addPoint(p5);
  box_side9.addPoint(p6);
  box_side9.addPoint(p7);
  box_side9.addPoint(p8);
  box_side9.setPixel();
  box_side9.setName("box_side9");

  //DrcSurfPolyFlat box_side10(box_side9);
  //box_side10.setName("box_side10");
  //Transform3D trans10(XYZVector(0,0,-50));
  //box_side10.addTransform(trans10);
  //box_side10.setReflectivity(refl_none);


  DrcOptVol box;
  box.setVerbosity(0);
  box.addSurface(box_side1);
  box.addSurface(box_side2);
  box.addSurface(box_side3);
  box.addSurface(box_side4);
  box.addSurface(box_side5);
  box.addSurface(box_side6);
  box.addSurface(box_side7);
  box.addSurface(box_side8);
  box.addSurface(box_side9);
  //box.addSurface(box_side10);
  box.addSurface(hole_side1);
  box.addSurface(hole_side2);
  box.addSurface(hole_side3);
  box.addSurface(hole_side4);
  box.addSurface(hole2_sphere2);
  box.addSurface(hole2_side1);
  box.addSurface(hole2_side2);
  box.addSurface(hole2_side3);
  box.addSurface(hole2_side4);
  box.setOptMaterial(vacuum);
  //box.setOptMaterial(quartz);
  box.setName("box");

  XYZPoint t5(-slab_width/2, +slab_height/2, -1000);
  XYZPoint t6(+slab_width/2, +slab_height/2, -1000);
  XYZPoint t7(-slab_width/2, -slab_height/2, -1000);
  XYZPoint t8(+slab_width/2, -slab_height/2, -1000);
  //s 5 6 at dist-1
  //  7 8
  //
  //t 5 6 at dist-1-1000
  //  7 8
  //
  DrcSurfPolyFlat slab_side1;
  slab_side1.setVerbosity(0);
  slab_side1.addPoint(s5); // face
  slab_side1.addPoint(s6);
  slab_side1.addPoint(s8);
  slab_side1.addPoint(s7);
  slab_side1.setName("slab_side1");
  DrcSurfPolyFlat slab_side2;
  slab_side2.addPoint(s5);//left
  slab_side2.addPoint(s7);
  slab_side2.addPoint(t7);
  slab_side2.addPoint(t5);
  slab_side2.setName("slab_side2");
  DrcSurfPolyFlat slab_side3;
  slab_side3.addPoint(s7);//down
  slab_side3.addPoint(s8);
  slab_side3.addPoint(t8);
  slab_side3.addPoint(t7);
  slab_side3.setName("slab_side3");
  DrcSurfPolyFlat slab_side4;
  slab_side4.addPoint(s6);//right
  slab_side4.addPoint(s8);
  slab_side4.addPoint(t8);
  slab_side4.addPoint(t6);
  slab_side4.setName("slab_side4");
  DrcSurfPolyFlat slab_side5;
  slab_side5.addPoint(s5);//top
  slab_side5.addPoint(s6);
  slab_side5.addPoint(t6);
  slab_side5.addPoint(t5);
  slab_side5.setName("slab_side5");
  DrcSurfPolyFlat slab_side6;
  slab_side6.addPoint(t5); // face
  slab_side6.addPoint(t6);
  slab_side6.addPoint(t8);
  slab_side6.addPoint(t7);
  slab_side6.setName("slab_side6");

  //slab_side1.setReflectivity(refl_perfect);
  slab_side2.setReflectivity(refl_perfect);
  slab_side3.setReflectivity(refl_perfect);
  slab_side4.setReflectivity(refl_perfect);
  slab_side5.setReflectivity(refl_perfect);
  slab_side6.setReflectivity(refl_perfect);

  DrcOptVol slab;
  slab.setVerbosity(0);
  slab.addSurface(slab_side1);
  slab.addSurface(slab_side2);
  slab.addSurface(slab_side3);
  slab.addSurface(slab_side4);
  slab.addSurface(slab_side5);
  slab.addSurface(slab_side6);
  slab.setOptMaterial(quartz);
  //slab.setOptMaterial(vacuum);
  slab.setName("slab");




  DrcOptDevSys opt_system;
  opt_system.addDevice(box);
  opt_system.addDevice(slab);
  opt_system.addDevice(lens);
  opt_system.addDevice(lens2);

  // couple surface 1 of device 1 with surface 2 of device 2
  //                       dev1   dev2    surf1     surf2
  opt_system.coupleDevice("lens","lens2","lens_sphere","lens2_sphere1");
  opt_system.coupleDevice("lens2","box","lens2_sphere2","hole2_sphere2");
  opt_system.coupleDevice("lens2","box","lens2_side1", "hole2_side1");
  opt_system.coupleDevice("lens2","box","lens2_side2", "hole2_side2");
  opt_system.coupleDevice("lens2","box","lens2_side3", "hole2_side3");
  opt_system.coupleDevice("lens2","box","lens2_side4", "hole2_side4");
  opt_system.coupleDevice("lens","box","lens_side1", "hole_side1");
  opt_system.coupleDevice("lens","box","lens_side2", "hole_side2");
  opt_system.coupleDevice("lens","box","lens_side3", "hole_side3");
  opt_system.coupleDevice("lens","box","lens_side4", "hole_side4");
  opt_system.coupleDevice("slab","lens","slab_side1", "lens_base");




  DrcOptDevManager* manager = new DrcOptDevManager();
  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(-450,-450,-300,450,450,500);"<<endl;
   geo<<"    view->SetRange(-450,-450,-450,450,450,450);"<<endl;
 geo<<"    Int_t i;"<<endl;
  geo<<"    view->SetView(180,45,90,i);"<<endl;
  geo<<"    view->SetView(90,0,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;
  
  // 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);
  

  //geo<<"}"<<endl;
  //exit(0);

  list<DrcPhoton> list_photon; // get list

  DrcPhoton ph;


  bool photons_exist = false;
  if (ioption==1)  
    {
      slab_width*=0.5;
      slab_height*=0.5;

      for (double xx=-slab_width/2; xx<=slab_width/2; xx+=slab_width/4)
	{
	  for (double yy=-slab_height/2; yy<=slab_height/2; yy+=slab_height/4)
	    {
	      ph.setPosition(XYZPoint(xx,yy,dist-500));
	      ph.setDirection(XYZVector(0,0,1));
	      ph.setWavelength(lambda);
	      list_photon.push_back(ph);
	    }
	}
      photons_exist = true;
      manager->setPhotonList(list_photon,"slab");
    }
  if (ioption==2)  
    {
      for (double angle1=-30; angle1<30.5; angle1+=5.0)
	{
	  for (double xx=-0.2*slab_width; xx<=0.21*slab_width; xx+=.40*slab_width)
	    {
	      double yy=0;
	      //for (double lambda1=330; lambda1<631; lambda1+=100)
	      {
		double angle;
		double lambda1 = lambda;
		ph.setPosition(XYZPoint(xx,yy,-500.5));
	      
		angle = pi/180.0 * angle1;
		ph.setDirection(XYZVector(sin(+angle),0,cos(+angle)));
		ph.setWavelength(lambda1);
		list_photon.push_back(ph);
	     
		   
	      }
	    }
	}
      photons_exist = true;
      manager->setPhotonList(list_photon,"slab");
    }
  if (ioption==3)
    {
      XYZPoint pos(0,-10,-45.0);
      XYZVector dir(0.0,1,-3); 
      double   beta = 0.80;
      photons_exist = manager->cerenkov(pos,dir,beta,100,200); // generate photons
    }
  if (ioption==4)  
    {
      TRandom ran;
      for (int ii=0; ii<30; ii++)
	{
	  ph.setPosition(XYZPoint(0,0,dist-20));
	  double phi =   ran.Uniform(0,2*pi);
	  double costh = ran.Uniform(-1,1);
	  double sinth = sqrt(1-costh*costh);
	  ph.setDirection(XYZVector(sinth*sin(phi),sinth*cos(phi),costh));
	  ph.setWavelength(lambda);
	  list_photon.push_back(ph);
	}
	
      photons_exist = true;
      manager->setPhotonList(list_photon,"slab");
    }
  if (ioption==5)
    {
      double angle = pi/180.0 * 0;
      ph.setDirection(XYZVector(0,sin(+angle),cos(+angle)));

      ph.setPosition(XYZPoint(0,5,-5));
      ph.setWavelength(400);
      list_photon.push_back(ph);
      ph.setPosition(XYZPoint(0,-5,-5));
      ph.setWavelength(400);
      list_photon.push_back(ph);


      ph.setPosition(XYZPoint(0,5,-5));
      ph.setWavelength(600);
      list_photon.push_back(ph);
      ph.setPosition(XYZPoint(0,-5,-5));
      ph.setWavelength(600);
      list_photon.push_back(ph);



      photons_exist = true;
      manager->setPhotonList(list_photon,"slab");
    }
  if (ioption==6)  
    {
      // straight lines.
      TRandom ran;
      for (double angle=5; angle<=30.5; angle+=5)
	{

	  cout<<" angle = "<<angle<<endl;
	  for (double lambda1=630; lambda1>329; lambda1-=10)
	    {
	      //double lambda1=lambda;
	      // go in x dir
	      double y1 = tan(angle*pi/180);
	      double x1 = y1;
	      double scale = 3/angle;
	      for (double y=-y1; y<= y1; y+= 2*y1/50*scale)
		{
		  double xx=ran.Uniform(-0.5*slab_width,0.5*slab_width);
		  double yy=ran.Uniform(-0.5*slab_height,0.5*slab_height);
		  ph.setPosition(XYZPoint(0,0,-0.5));
		  double x = x1;
		  ph.setDirection(XYZVector(x,y,1).Unit());
		  ph.setWavelength(lambda1);
		  list_photon.push_back(ph);
		  ph.setDirection(XYZVector(-x,y,1).Unit());
		  ph.setWavelength(lambda1);
		  list_photon.push_back(ph);
		}
	      for (double x=-x1; x<= x1; x+= 2*y1/50*scale)
		{
		  double xx=ran.Uniform(-0.5*slab_width,0.5*slab_width);
		  double yy=ran.Uniform(-0.5*slab_height,0.5*slab_height);
		  ph.setPosition(XYZPoint(0,0,-0.5));
		  double y = y1;
		  ph.setDirection(XYZVector(x,y,1).Unit());
		  ph.setWavelength(lambda1);
		  list_photon.push_back(ph);
		  ph.setDirection(XYZVector(x,-y,1).Unit());
		  ph.setWavelength(lambda1);
		  list_photon.push_back(ph);
		}
	     
	    }
	}
      photons_exist = true;
      manager->setPhotonList(list_photon,"slab");
    }

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

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




  geo<<"}"<<endl;
  //
  // the intention is to play around with routines.

  // create a list of photons in bar
  //exit(1);
  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]\",600,-300,300);"<<endl;
  scr<<"    hgr->SetMarkerStyle(7);"<<endl;
  scr<<"    hgr->SetMarkerSize(0.5);"<<endl;
  scr<<"    hgr->SetMinimum(-300);"<<endl;
  scr<<"    hgr->SetMaximum(+300);"<<endl;
  scr<<"    hgr->Draw(\"POL\");"<<endl;
  int icnt_measured = 0;
  int icnt_flying   = 0;
  int icnt_lost     = 0;
  int icnt_absorbed = 0;

  //TRotation rotInv_screen = rot_screen.Inverse();
  list<DrcPhoton>::iterator iph;
  DrcPhoton ph_old;

  double dista,distb;
  for(iph=list_photon.begin(); iph != list_photon.end(); ++iph) 
    {


      //(*iph).print();
      if ((*iph).fate()==Drc::PhotMeasured)
	{
	  icnt_measured++;
	  XYZPoint pos = (*iph).position();
	  //pos -= shift_screen;
	  //pos.Transform(rotInv_screen);

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

	  //cout<<" i,p,d="<<icnt_measured<<" "
	  //<<(*iph).position()<<" "<<(*iph).direction()<<endl;

	  if ((icnt_measured%2)==0)
	    {
	      //cout<<icnt_measured<<endl;
	      XYZPoint p1((*iph).position());
	      XYZVector d1((*iph).direction());
	      XYZPoint p2(ph_old.position());
	      XYZVector d2(ph_old.direction());
	      double d1x = d1.X();
	      double d1z = d1.Z();
	      double d2x = d2.X();
	      double d2z = d2.Z();
	      double p1x = p1.X();
	      double p1z = p1.Z();
	      double p2x = p2.X();
	      double p2z = p2.Z();


	      double d = (d1x - d1z/d2z*d2x);
	      double lam = (p2x-p1x)/d;
	      //cout<<p1x<<" "<<p1z<<" "<<d1x<<" "<<d1z<<endl;
	      //cout<<p2x<<" "<<p2z<<" "<<d2x<<" "<<d2z<<endl;
	      
	      //cout<<" "<<(p2x-p1x)<<" "<<d<<" "<<lam*d1z<<endl;
	      cout<<" distz = "<<p1z+lam*d1z<<endl;
	      //cout<<" distx = "<<p1x+lam*d1x<<endl;

	      if (icnt_measured==2) dista = p1z+lam*d1z;
	      if (icnt_measured==14) distb = p1z+lam*d1z;

	    }

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

}