#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 <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 "PndDrcPhoton.h"
#include "PndDrcOptReflAbs.h"
#include "PndDrcSurfAbs.h"
#include "PndDrcOptMatAbs.h"
#include "PndDrcOptDev.h"
#include "PndDrcOptDevSys.h"
#include "PndDrcOptDevManager.h"
#include "PndDrcUtil.h"



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

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

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

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

  map<string,PndDrcOptDev*>::iterator idev = m_mapDev.find(str);
  if (idev != m_mapDev.end())
    {
      dev = (*idev).second;
    }
  else
    {
	  cerr<<" *** PndDrcOptDevManager::device: device not in map: "
	      <<str<<endl;
	  exit(EXIT_FAILURE);
    }
  return dev;
}
//----------------------------------------------------------------------
 void PndDrcOptDevManager::addDeviceSystem(const PndDrcOptDevSys& sys)
{
  PndDrcOptDevSys* tmp = new PndDrcOptDevSys();
  *tmp = sys;
  m_listDevSys.push_back(tmp);


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

  list<PndDrcOptDev*>::const_iterator idev;
  list<string>::const_iterator        isys;
  list<int>::const_iterator           isys_copy;
  
  
  cout<<" size="<<tmp->listSysOri().size()<<endl;
  
  int ii=0;
   for(list<string>::const_iterator isys1  = tmp->listSysOri().begin();
       isys1 != tmp->listSysOri().end();
       ++isys1) 
    {
      string str = *isys1;
      
      cout<<++ii<<" "<<str<<endl;
      if (ii==10) exit(1);
      
      
    }
   

  ii=0;
  
  for(idev        = tmp->deviceList().begin(),
	isys      = tmp->listSysOri().begin(),
	isys_copy = tmp->listSysOriCopy().begin();
      idev != tmp->deviceList().end(); 
      ++idev,
	++isys,
	++isys_copy) 
    {
      cout<<(++ii)<<" "<<(*idev)->name()<<endl;
      //cout<<(++ii)<<" "<<icopy<<endl;
      cout<<(ii)<<" "<<(*isys)<<endl;
      //cout<<(++ii)<<" "<<(isys_copy)<<endl;
      
      string str1 = (*idev)->name() + "_" + itoa((*idev)->copyNumber(),6);
      //string str2 = (*isys);// + "_" + itoa((*isys_copy),6);
      string str  =  str1;//+str2;
      
      pair<string,PndDrcOptDev*> p(str,(*idev));
      pair< map<string,PndDrcOptDev*>::iterator, bool > idev1 = m_mapDev.insert(p);

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

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


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


  if (d1 && d2)
    {
      list<PndDrcSurfAbs*>::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<<" *** PndDrcOptDevManager::coupleDevice: surface not found"<<endl;
	  if (!s1) cerr<<"     "<<surf1<<" missing"<<endl;
	  if (!s2) cerr<<"     "<<surf2<<" missing"<<endl;
	  exit(EXIT_FAILURE);
	}
    }
  else
    {
      cerr<<" *** PndDrcOptDevManager::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 PndDrcOptDevManager::print(fstream& stream) const
{
  list<PndDrcOptDevSys*>::const_iterator isys;

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

  HitMapType    hit_map;
  VolMapType    vol_map;
  SysOriMapType sys_ori_map;
  HitMapType::iterator    hit_map_iter;
  VolMapType::iterator    vol_map_iter;
  SysOriMapType::iterator sys_ori_map_iter;

  list<PndDrcOptDevSys*>::const_iterator isys;
  list<PndDrcOptDev*>::const_iterator    idev;
  list<PndDrcSurfAbs*>::const_iterator   isurf;

  list<string>::const_iterator           isys_ori;
  list<int>::const_iterator              isys_ori_copy;
  
  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)
	{
	  PndDrcOptDevSys* sys = (*isys);
	  list<PndDrcOptDev*> list_dev = sys->deviceList();
	  
	  for (idev=list_dev.begin(),
		 isys_ori=sys->listSysOri().begin(),
		 isys_ori_copy=sys->listSysOriCopy().begin(); 
	       idev != list_dev.end(); 
	       idev++,
		 isys_ori++,
		 isys_ori_copy++)
	    {
	      PndDrcOptDev* dev = (*idev);
	      list<PndDrcSurfAbs*> list_surf = dev->surfaceList();
	      XYZPoint pos_new;

	      PndDrcSurfAbs* surf_closest=0;
	      if (dev->radiator())
		{
		  for (isurf=list_surf.begin(); isurf != list_surf.end(); isurf++)
		    {
		      PndDrcSurfAbs* 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)
		{
		  PndDrcSurfAbs* 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) );
			  pair<string,int> sys_ori_pair;
			  sys_ori_pair.first  = (*isys_ori);
			  sys_ori_pair.second = (*isys_ori_copy);
			  sys_ori_map.insert(SysOriValuePair(path_length,sys_ori_pair));
			}	      
		    }
		}
	    }
	} 
    } // while hit







  if (vol_map.size()>0)
    {

      // map is already sorted.

      p_in.clear();
      p_out.clear();
      vol_name.clear();
      vol_copy.clear();
      sys_name.clear();
      sys_copy.clear();

      for (vol_map_iter=vol_map.begin(),
	     hit_map_iter=hit_map.begin(),
	     sys_ori_map_iter=sys_ori_map.begin(); 
	   vol_map_iter != vol_map.end(); 
	   ++ vol_map_iter, 
	     ++ hit_map_iter,
	     ++ sys_ori_map_iter)
	{
	  pair<XYZPoint,XYZPoint> r = (*hit_map_iter).second;
	  pair<string,int> v = (*vol_map_iter).second;
	  pair<string,int> s = (*sys_ori_map_iter).second;
	  
	  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);
	      vol_copy.push_back(v.second);
	      sys_name.push_back(s.first);
	      sys_copy.push_back(s.second);
	    }
	}
      
      return (p_in.size() >0) ? true : false;
    }
  
  return false;
}
//----------------------------------------------------------------------
bool PndDrcOptDevManager::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>      vol_copy_list;
  list<string>   sys_name_list;
  list<int>      sys_copy_list;

  
  if(intersect(pos,dir,
	       p_in_list,p_out_list,
	       vol_name_list,sys_name_list,
	       vol_copy_list,sys_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_vol_copy;
      list<string>::const_iterator   it_sys_name;
      list<int>::const_iterator      it_sys_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_vol_copy = vol_copy_list.begin(),
	     it_sys_name = sys_name_list.begin(),
	     it_sys_copy = sys_copy_list.begin();
	   it_p_in != p_in_list.end();
	   ++it_p_in,
	     ++it_p_out,
	     ++it_vol_name,
	     ++it_vol_copy,
	     ++it_sys_name,
	     ++it_sys_copy)
	{
	  XYZPoint     r1 = (*it_p_in);
	  XYZPoint     r2 = (*it_p_out);
	  string vol_name = (*it_vol_name);
	  int    vol_copy = (*it_vol_copy);
	  string sys_name = (*it_sys_name);
	  int    sys_copy = (*it_sys_copy);
	  

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

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

	}
      
      return result;
    }

  return true;
}
//----------------------------------------------------------------------
bool PndDrcOptDevManager::cerenkov(const string& vol_name, const string& sys_name,
				   int vol_copy, int sys_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;
      if (lam1==lam2) lam2=lam1+1;
      nu_mean = (log(lam2)-log(lam1))/(1.0/lam1 - 1.0/lam2);
      if (m_verbosity>=3) cout<<" PndDrcOptDevManager::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<<"*   : PndDrcOptDevManager::cerenkov: beta="<<beta<<endl;
      cerr<<"      reset beta=0.9999"<<endl;
      beta = 0.9999;
    }


  PndDrcOptDev* dev=device(vol_name,sys_name,vol_copy,sys_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);


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

    }
     
  return true;

}
//----------------------------------------------------------------------
void PndDrcOptDevManager::setPhotonList(list<PndDrcPhoton>& photon_list,
					string vol_name, string sys_name, 
					int ivol_copy, int isys_copy)
{
  m_listPhoton = photon_list;

  // find the internal device pointer
  PndDrcOptDev* dev=device(vol_name,sys_name,ivol_copy,isys_copy);
  
  list<PndDrcPhoton>::iterator iph;
  for (iph = m_listPhoton.begin(); iph != m_listPhoton.end(); ++iph)
  {
    (*iph).setDevice(dev);
  }
}
//----------------------------------------------------------------------
void PndDrcOptDevManager::propagate()
{
  list<PndDrcPhoton>::iterator iph;

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

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

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

  list<PndDrcOptDevSys*>::const_iterator isys;


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