#include <iostream>
using std::cout;
using std::cerr;
using std::cin;
using std::endl;
using std::hex;

#include <valarray>
using std::valarray;

#include <vector>
using std::vector;

#include <string>
using std::string;

#include <list>
using std::list;

#include <map>
using std::map;

#include <limits>
using std::numeric_limits;

#include <fstream>
using std::fstream;

#include <utility>
using std::pair;

//#include <cmath>

#include "TObject.h"
#include "TVector3.h"
#include "TRandom.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 "DrcOptMatAbs.h"
#include "DrcOptDev.h"
#include "DrcOptDevSys.h"
#include "DrcOptDevManager.h"
#include "DrcUtil.h"



typedef map< double, pair<XYZPoint,XYZPoint> > HitMapType;
typedef map< double, pair<string,int> >        VolMapType;
typedef HitMapType::value_type HitValuePair;
typedef VolMapType::value_type VolValuePair;

// pointer to the object made by Construct()
DrcOptDevManager* DrcOptDevManager::m_instance = 0;

//----------------------------------------------------------------------
DrcOptDevManager::DrcOptDevManager()
{
  if(m_instance)
    {
      cerr<<" There is only one DrcOptDevManager allowed! "<<endl;
      exit(EXIT_FAILURE);
    }
  m_verbosity = 0;
}
//----------------------------------------------------------------------
DrcOptDevManager::~DrcOptDevManager()
{
  list<DrcOptDevSys*>::const_iterator isys;

  for(isys=m_listDevSys.begin(); isys != m_listDevSys.end(); ++isys) 
    {
      delete (*isys);
   }
}
//----------------------------------------------------------------------
DrcOptDev* DrcOptDevManager::device(string s,int icopy)
{
  string     str = s + "_" + itoa(icopy,6);
  DrcOptDev* dev = 0;

  map<string,DrcOptDev*>::iterator idev = m_mapDev.find(str);
  if (idev != m_mapDev.end())
    {
      dev = (*idev).second;
    }
  else
    {
	  cerr<<" *** DrcOptDevManager::device: device not in map: "
	      <<str<<endl;
	  exit(EXIT_FAILURE);
    }
  return dev;
}
//----------------------------------------------------------------------
 void DrcOptDevManager::addDeviceSystem(const DrcOptDevSys& sys, int icopy)
{

  if (icopy>99999)
    {
	  cerr<<" *** DrcOptDevManager::addDeviceSystem: sorry, copy limited to"
	      <<" 5 digits."<<endl;
	  exit(EXIT_FAILURE);
    }


  DrcOptDevSys* tmp = new DrcOptDevSys();
  *tmp = sys;
  tmp->setCopyNumber(icopy);
  m_listDevSys.push_back(tmp);


  // add for fast finding of devices per name the combined name and copy number to
  // a map;

  list<DrcOptDev*>::iterator idev;
  for(idev  = tmp->deviceList().begin(); idev != tmp->deviceList().end(); ++idev) 
    {
      string str = (*idev)->name() + "_" + itoa(icopy,6);

      pair<string,DrcOptDev*> p(str,(*idev));
      pair<map<string,DrcOptDev*>::iterator,bool> idev1 = m_mapDev.insert(p);

      if (!idev1.second)
	{
	  cerr<<" *** DrcOptDevManager::addDeviceSystem: device already present \n"
	      <<"     name,copy: "<<tmp->name()<<" "<<icopy<<endl;

	  exit(EXIT_FAILURE);
	}
      else
	{
	  //cout<<" inserted "<<str<<endl;
	}
    }
  


}
//----------------------------------------------------------------------
void DrcOptDevManager::coupleDevice(string dev1,    string dev2,
				    string surf1,   string surf2,
				    int    copy1,   int    copy2)
{
  list<DrcOptDevSys*>::const_iterator isys;
  DrcOptDev*  d1=device(dev1,copy1);
  DrcOptDev*  d2=device(dev2,copy2);
  DrcSurfAbs* s1=0;
  DrcSurfAbs* s2=0;


  if (d1 && d2)
    {
      list<DrcSurfAbs*>::const_iterator isurf;
      for(isurf=(d1->surfaceList()).begin(); 
	  isurf != (d1->surfaceList()).end(); ++isurf) 
	{
	  if ((*isurf)->name()== surf1) s1=(*isurf); 
	}
      for(isurf=d2->surfaceList().begin(); isurf != d2->surfaceList().end(); ++isurf) 
	{
	  if ((*isurf)->name()== surf2) s2=(*isurf); 
	}
      if (s1 && s2)
	{
	  s1->setCoupled(d2,s2);
	  s2->setCoupled(d1,s1);
	}
      else
	{
	  cerr<<" *** DrcOptDevManager::coupleDevice: surface not found"<<endl;
	  if (!s1) cerr<<"     "<<surf1<<" missing"<<endl;
	  if (!s2) cerr<<"     "<<surf2<<" missing"<<endl;
	  exit(EXIT_FAILURE);
	}
    }
  else
    {
      cerr<<" *** DrcOptDevManager::coupleDevice: device not found"<<endl;
      if (!d1) cerr<<"     "<<dev1<<" missing (copy="<<copy1<<")"<<endl;
      if (!d2) cerr<<"     "<<dev2<<" missing (copy="<<copy2<<")"<<endl;
      exit(EXIT_FAILURE);
    }
}
//----------------------------------------------------------------------
void DrcOptDevManager::print(fstream& stream) const
{
  list<DrcOptDevSys*>::const_iterator isys;

  for(isys=m_listDevSys.begin(); isys != m_listDevSys.end(); ++isys)
    {
      (*isys)->print(stream);
    } 
}
//----------------------------------------------------------------------
bool DrcOptDevManager::intersect(const XYZPoint&  pos, 
				 const XYZVector& dir,
				 list<XYZPoint>& p_in, 
				 list<XYZPoint>& p_out, 
				 list<string>&    vol_name,
				 list<int>&       copy,
				 double           range) const
{
  const double eps = 1.0e-9;
  DrcPhoton ph;
  ph.setPosition(pos);
  ph.setDirection(dir.Unit());

  HitMapType hit_map;
  VolMapType vol_map;
  HitMapType::iterator hit_map_iter;
  VolMapType::iterator vol_map_iter;

  list<DrcOptDevSys*>::const_iterator isys;
  list<DrcOptDev*>::const_iterator    idev;
  list<DrcSurfAbs*>::const_iterator   isurf;

  double path_length;
  bool hit;
  
  pair<XYZPoint,XYZPoint> hit_pair;
  bool hit1=true;
  int ihit = 0;
  while (hit1 && (ph.position()-pos).Mag2()<range*range )
    {
      hit1 = false;
	      double path_shortest = numeric_limits<double>::max();
      for (isys=m_listDevSys.begin(); isys != m_listDevSys.end(); ++isys)
	{
	  DrcOptDevSys* sys = (*isys);
	  list<DrcOptDev*> list_dev = sys->deviceList();

	  for (idev=list_dev.begin(); idev != list_dev.end(); idev++)
	    {
	      DrcOptDev* dev = (*idev);
	      list<DrcSurfAbs*> list_surf = dev->surfaceList();
	      XYZPoint pos_new;

	      DrcSurfAbs* surf_closest=0;
	      if (dev->radiator())
		{
		  for (isurf=list_surf.begin(); isurf != list_surf.end(); isurf++)
		    {
		      DrcSurfAbs* surf = (*isurf);
		      hit = surf->surfaceHit(ph,pos_new,path_length);
		      if (hit)
			{
			  if (path_length < path_shortest)
			    {
			      path_shortest = path_length;
			      surf_closest  = surf;
			    }
			}
		    }
		}

	      if (surf_closest)
		{
		  DrcSurfAbs* surf = surf_closest;
		  hit = surf->surfaceHit(ph,pos_new,path_length);
		  if (hit) hit1=true;

		  if (hit)
		    {
		      ph.setPosition(pos_new+ph.direction()*2*eps); // bring it in.
		      ihit++;

		      //cout<<" ihit%2="<<ihit%2<<endl;
		      //cout<<pos_new.X()<<" "<<pos_new.Y()<<" "<<pos_new.Z()<<endl;//###

		      if (ihit%2 == 1) 
			{
			  hit_pair.first  = pos_new;
			}
		      else if (ihit%2 == 0) 
			{
			  hit_pair.second = pos_new;
	
			  hit_map.insert( HitValuePair(path_length,hit_pair) );
			  pair<string,int> vol_pair;
			  vol_pair.first = dev->name();
			  vol_pair.second = dev->copyNumber();
			  vol_map.insert( VolValuePair(path_length,vol_pair) );
			}	      
		    }
		}
	    }
	} 
    } // while hit







  if (vol_map.size()>0)
    {

      // map is already sorted.

      p_in.clear();
      p_out.clear();
      vol_name.clear();
      copy.clear();

      for (vol_map_iter=vol_map.begin(),
	     hit_map_iter=hit_map.begin(); 
	   vol_map_iter != vol_map.end(); 
	   ++vol_map_iter, ++ hit_map_iter)
	{
	  pair<XYZPoint,XYZPoint> r = (*hit_map_iter).second;
	  pair<string,int> v = (*vol_map_iter).second;
	  //double path_length = (*vol_map_iter).first;
	  //cout<<" path_length = "<<path_length<<endl;
	  if ((r.first-r.second).Mag2() > 1.0e-18) // no flat devices.
	    {
	      p_in.push_back(r.first);
	      p_out.push_back(r.second);
	      vol_name.push_back(v.first);
	      copy.push_back(v.second);
	    }
	}
      
      return (p_in.size() >0) ? true : false;
    }
  
  return false;
}
//----------------------------------------------------------------------
bool DrcOptDevManager::cerenkov(const XYZPoint& pos, const XYZVector& dir, 
				double beta, int n, double range, 
				double nu1, double nu2)
{

  list<XYZPoint> p_in_list;
  list<XYZPoint> p_out_list;
  list<string>   vol_name_list;
  list<int>      copy_list;

  
  if(intersect(pos,dir,p_in_list,p_out_list,vol_name_list,copy_list,range))
    {
      list<XYZPoint>::const_iterator it_p_in;
      list<XYZPoint>::const_iterator it_p_out;
      list<string>::const_iterator   it_vol_name;
      list<int>::const_iterator      it_copy;
	  
      bool result = false;

      for (it_p_in = p_in_list.begin(),
	     it_p_out = p_out_list.begin(),
	     it_vol_name = vol_name_list.begin(),
	     it_copy = copy_list.begin();
	   it_p_in != p_in_list.end();
	   ++it_p_in,
	     ++it_p_out,
	     ++it_vol_name,
	     ++it_copy)
	{
	  XYZPoint     r1 = (*it_p_in);
	  XYZPoint     r2 = (*it_p_out);
	  string vol_name = (*it_vol_name);
	  int copy_number = (*it_copy);


	  //cout<<" r1= "<<r1.X()<<" "<<r1.Y()<<" "<<r1.Z()<<endl;//###
	  //cout<<" r2= "<<r2.X()<<" "<<r2.Y()<<" "<<r2.Z()<<endl;//###

	  if (cerenkov(vol_name,copy_number,r1,r2,beta,n,nu1,nu2)) result = true;
	  

	}
      
      return result;
    }

  return true;
}
//----------------------------------------------------------------------
bool DrcOptDevManager::cerenkov(const string& vol_name, int copy, 
				const XYZPoint& r1, const XYZPoint& r2, 
				double beta, int inum, double lam1, double lam2)
{
  static bool first=true;

  //const double cc = 29.9792;      // c in cm/ns
  const double pi = 3.1415926535;

  static double nu_mean,fn0;


  if (first)
    {
      first = false;
      nu_mean = (log(lam2)-log(lam1))/(1.0/lam1 - 1.0/lam2);
      if (m_verbosity>=3) cout<<" DrcOptDevManager::cerenkov nu_mean = "
			      <<nu_mean<<endl;

      double kk   = 370;           // cm-1 eV-1, PRD 50 (1994) 1261
      double c    = 2.99792458e10; // cm/s
      double hbar = 6.5821220e-16; // eV s
      
      double nu1  = c/(lam1*1.0e-7); // s-1 lower end of visible spectrum
      double nu2  = c/(lam2*1.0e-7); // s-1 upper end
  
      fn0 = kk*(hbar*2*pi)*(nu1-nu2); // cm-1
      fn0 *= 0.1;                     // mm-1 
      if (m_verbosity>=3) 
	{
	  cout<<" N0 = "<<fn0<<" mm-1"
	      <<" for nu1,nu2 = "<<lam1<<"-"<<lam2<<endl;
	  cout<<"                that is dE = "<<(hbar*2*pi)*(nu1-nu2)<<" eV"<<endl;
	}
    }



  if (beta>1.001) 
    {
      cerr<<"*   : DrcOptDevManager::cerenkov: beta="<<beta<<endl;
      cerr<<"      reset beta=0.9999"<<endl;
      beta = 0.9999;
    }


  DrcOptDev* dev=device(vol_name,copy);


  // get number of photons to sample
  double costh = 1.0/((dev->optMaterial().refIndex(nu_mean))*beta);
  if (costh>1) return 0; // particle too slow
  

  //cout<<"inum="<<fn0<<" "<<(r1-r2).Mag()<<" "<<(1-costh*costh)<<endl;//###


  int inum1 = inum;
  if (inum==0) // 0 = Cherenkov distributed
    {
      double path = sqrt((r1-r2).Mag2());
      inum1 = int(fn0*path*(1-costh*costh) + 0.5);
    }


  if (inum1==0) return false; // no photons.

  double   x1,x2,x,lambda;
  XYZVector diff;
  XYZPoint  r;
  XYZVector beta_ph;
  double    dtheta,dphi;


  //cout<<" generate photons: "<<inum1<<endl;

  for (int i=0; i<inum1; i++) // loop over photons
    {
      do
	{
	  // sample a lambda distribution with a 1/lambda^2 shape. 
	  x1     = 1.0/lam2;               // inv. lambda start at x1
	  x2     = 1.0/lam1;               // inv. lambda end at x2
	  x      = m_ran.Uniform(x1,x2);  // inv. lambda random range
	  lambda = 1.0/x;                        // yields 1/lambda^2 dist.
	  //cout<<" x,lambda="<<x<<" "<<lambda<<endl;//###
	  double n = dev->optMaterial().refIndex(lambda);
	  costh = 1.0L/(beta*n);
	} while (costh>1); 

      diff = r2 - r1;
      r    = r1 + m_ran.Uniform(1.0)*diff;
      beta_ph   = diff.Unit();      // velocity vector with length 1
      dphi      = m_ran.Uniform(0,2*pi); 


      //cout<<dphi<<endl;//###


      dtheta = acos(costh);

      beta_ph   = angle_tra(beta_ph,dtheta,dphi);


      DrcPhoton ph;
      ph.setPosition(r);
      ph.setDirection(beta_ph);
      ph.setWavelength(lambda);
      ph.setDevice(dev);
      m_listPhoton.push_back(ph);

    }
     
  return true;

}
//----------------------------------------------------------------------
void DrcOptDevManager::setPhotonList(list<DrcPhoton>& photon_list,
				     string vol_name, int copy)
{
  m_listPhoton = photon_list;

  // find the internal device pointer
  DrcOptDev* dev=device(vol_name,copy);
  
  list<DrcPhoton>::iterator iph;
  for (iph = m_listPhoton.begin(); iph != m_listPhoton.end(); ++iph)
  {
    (*iph).setDevice(dev);
  }
}
//----------------------------------------------------------------------
void DrcOptDevManager::propagate()
{
  list<DrcPhoton>::iterator iph;

  for (iph = m_listPhoton.begin(); iph != m_listPhoton.end(); ++iph)
    {
      DrcOptDev* dev = (*iph).device();
      dev->propagate((*iph));
    }
}
//----------------------------------------------------------------------
void DrcOptDevManager::addTransform(const Transform3D& trans)
{
  list<DrcOptDevSys*>::const_iterator isys;

  for (isys=m_listDevSys.begin(); isys != m_listDevSys.end(); ++isys)
    {
      DrcOptDevSys* sys = (*isys);
      sys->addTransform(trans);
    }
}
//----------------------------------------------------------------------
void DrcOptDevManager::print()
{

  cout<<"---------------------------- Start DrcOptDevManager objects and pointers"
      <<endl;

  list<DrcOptDevSys*>::const_iterator isys;


  for (isys=m_listDevSys.begin(); isys != m_listDevSys.end(); ++isys)
    {
      cout<<"\n\n"<<hex<<int(*isys)<<" "<<(*isys)->name()<<endl;
      (*isys)->print();
    }
  cout<<"---------------------------- End  DrcOptDevManager objects and pointers"
      <<endl;
}