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

#include <valarray>
using std::valarray;

#include <fstream>
using std::fstream;

#include <string>
using std::string;

#include <list>
using std::list;

#include "TObject.h"
#include "TVector3.h"
#include "TRandom.h"

#include "DrcOptMatAbs.h"
#include "DrcOptMatLithotecQ0.h"


//----------------------------------------------------------------------
DrcOptMatLithotecQ0::DrcOptMatLithotecQ0()
{
  m_f1 = 28.17847;
  m_e1 = 10.28086;
  m_f2 = 207.5092;
  m_e2 = 15.77682;
}
//----------------------------------------------------------------------
DrcOptMatLithotecQ0* DrcOptMatLithotecQ0::clone() const
{
  return new DrcOptMatLithotecQ0(*this);
}
//----------------------------------------------------------------------
void DrcOptMatLithotecQ0::copy(const DrcOptMatLithotecQ0& mat)
{
  m_f1  = mat.m_f1;
  m_e1  = mat.m_e1;
  m_f2  = mat.m_f2;
  m_e2  = mat.m_e2;
  m_ran = mat.m_ran;
}//----------------------------------------------------------------------
DrcOptMatLithotecQ0::DrcOptMatLithotecQ0(const DrcOptMatLithotecQ0& mat)
  : DrcOptMatAbs(mat)
{
  if (mat.m_verbosity>=1) cout<<"  DrcOptMatLithotecQ0::DrcOptMatLithotecQ0"
			    <<"(const DrcOptMatLithotecQ0&) "  
			    <<mat.m_name<<endl;
  copy(mat);
} 
//----------------------------------------------------------------------
DrcOptMatLithotecQ0& DrcOptMatLithotecQ0::operator=(const DrcOptMatLithotecQ0& mat)
{
  if (mat.m_verbosity>=1) cout<<"  DrcOptMatLithotecQ0::operator="
			    <<"(const DrcOptMatLithotecQ0&) "  
			    <<mat.m_name<<endl;
  if (&mat != this)
    {
      static_cast<DrcOptMatAbs&>((*this)) = mat; // assignment of base class part.
      copy(mat);
    }
  return *this;
}
//----------------------------------------------------------------------
double DrcOptMatLithotecQ0::refIndex(const double lambda) const
{
  static const double pi    = 3.1415926535;
  static const double hbarc = 197.3269602;  // Mev fm

  if (lambda<0) return 1.47; // average value.

  //cout<<" DrcOptMatLithotecQ0::refIndex: lambda = "<<lambda<<endl;//###
  double e = 2*pi*hbarc/lambda; // energy [eV]
  
  if (e<m_e1)
    {
      // refraction index
      return sqrt(1.0L + m_f1/(m_e1*m_e1-e*e) + m_f2/(m_e2*m_e2-e*e));
    }     
  else
    {
      cerr<<" *** DrcOptMatLithotecQ0::refIndex: too short wavelength. "
	  <<" lambda = "<<lambda<<" ... abort"<<endl;
      exit(EXIT_FAILURE);
    }

  return 1.47; 
}  
//----------------------------------------------------------------------
double DrcOptMatLithotecQ0::refIndexDeriv(const double lambda) const
{
  static const double pi    = 3.1415926535;
  static const double hbarc = 197.3269602;  // Mev fm

  double e = 2*pi*hbarc/lambda; // energy [eV]
  
  if (e<m_e1)
    {
      // refraction index
      return e/refIndex(lambda) * ( m_f1/(m_e1*m_e1-e*e)/(m_e1*m_e1-e*e) + 
				    m_f2/(m_e2*m_e2-e*e)/(m_e2*m_e2-e*e)  ) 
	* (-2*pi*hbarc/lambda/lambda);
    }     
  else
    {
      cerr<<" *** DrcOptMatLithotecQ0::refIndexDeriv: too short wavelength. "
	  <<" lambda = "<<lambda<<" ... abort"<<endl;
      exit(EXIT_FAILURE);
    }

 return 0;
}
//----------------------------------------------------------------------
bool DrcOptMatLithotecQ0::absorptionFlag(double lambda, double length) const
{
  // Rayleigh scattering. 
  static const double clarity = 2100*1000; // @ 633 nm in mm (2100 m) (=278m at 400nm)

  double trans = exp(-(length)/(clarity*pow(lambda/663,4)));
  double cmp   = m_ran.Uniform(1.0);


  if (cmp>trans)
    {
      return true; // absorbed
    } 

  return false; // no absorption.
}