//
// $Id: Material.cpp,v 1.6 2001/01/25 14:22:18 hobbs Exp $
//
// File: material.cpp
// Purpose: Represent a material
// Created: 03-SEP-1998   Harry L. Melanson
//
// $Revision: 1.6 $
//
//

////////////////////////////////////////////////////////////////////////
//
// Include files
//
#include <cmath>
#include <iomanip>
#include "material/Material.hpp"

using namespace std;

namespace D0Material {

  ////////////////////////////////////////////////////////////////////////
  //
  // Constructor to generate a material, calculating radiation length using (A,Z)
  //
  Material::Material(const string& name, 
		     const string& symbol,
		     double zeff, 
		     double aeff,
		     double rho) :
    d_name(name), 
    d_symbol(symbol),
    d_zeff(zeff),
    d_aeff(aeff),
    d_rho(rho)
    {
      ComputeLradTsaiFactor();
    }
  ////////////////////////////////////////////////////////////////////////
  //
  // Constructor to generate a material, given (A,Z) and radiation length
  //
  Material::Material(const string& name, 
		     const string& symbol,
		     double zeff, 
		     double aeff,
		     double rho,
		     double x0) :
    d_name(name), 
    d_symbol(symbol),
    d_zeff(zeff),
    d_aeff(aeff),
    d_rho(rho),
    d_x0(x0)
    {}
  
  ////////////////////////////////////////////////////////////////////////
  //
  void Material::ComputeLradTsaiFactor()
    {
      //
      //  Compute Tsai's Expression for the Radiation Length
      //  (Phys Rev. D50 3-1 (1994) page 1254)
      //
      //  1/X_0 = 4 alpha r_e^2 N_A/A { Z^2 {L_RAD - f(Z)] + Z L_RAD'}
      //
      //  Fills data member d_X0, in units of gm/cm**2
      //
      
      //
      //  First calculate f(z)
      //
      const double fine_structure_constant = 1.0 / 137.0359895;
      
      const double a = fine_structure_constant * d_zeff;
      const double a2 = a * a;
      const double a4 = a2 * a2;
      
      const double fz = a2 * ((1.0 / (1.0 + a2))
			      + 0.20206
			      - 0.03690 * a2
			      + 0.00830 * a4
			      - 0.00200 * a2* a4);
      //
      //  Then calculate L_RAD and L_RAD'
      //
      double Lrad, Lradp;
      
      const double Lrad_light[]  = {5.31  , 4.79  , 4.74 ,  4.71} ;
      const double Lradp_light[] = {6.144 , 5.621 , 5.805 , 5.924} ;
      
      double logZ3 = 0.0;
      if( d_zeff>0.0 ) logZ3=log(d_zeff)/3.;
      
      const int iz = (int)(d_zeff+0.5) - 1 ;
      
      if (iz <= 3) {                           //  H, He, Li, Be:
	Lrad = Lrad_light[iz];                 //  Look up values
	Lradp = Lradp_light[iz]; 
      }
      else {                                   //  Heavier:
	Lrad = log(184.15) - logZ3;            //  Calculate values
	Lradp = log(1194.) - 2*logZ3; 
      }
      //
      //  Finally, put it all together
      //
      const double factor = 1.0 / 716.408;     // factor == (4 alpha r_e^2 N_A)

      if (d_aeff != 0.0)
	d_x0 = (factor / d_aeff) * d_zeff * (d_zeff * (Lrad - fz) + Lradp);
      else
	d_x0 = 0.0;
      
      //
      //  Convert to X_0 (in units of gm/cm**2)
      //
      if (d_x0 != 0.0)
	d_x0 = 1.0 / d_x0;
      
    }
  
  
  ////////////////////////////////////////////////////////////////////////
  //
  //  Copy constructor
  //
  Material::Material(const Material &right)
    {      
      d_name   = right.d_name;
      d_symbol = right.d_symbol;
      d_aeff   = right.d_aeff;
      d_zeff   = right.d_zeff;
      d_rho    = right.d_rho;
      d_x0     = right.d_x0;
    }
  
  
  ////////////////////////////////////////////////////////////////////////
  //
  const Material& Material::operator=(const Material &right)
    {if (&right != this) {
      d_name   = right.d_name;
      d_symbol = right.d_symbol;
      d_aeff   = right.d_aeff;
      d_zeff   = right.d_zeff;
      d_rho    = right.d_rho;
      d_x0     = right.d_x0;
	}
    return *this;
    }
  
  
  ////////////////////////////////////////////////////////////////////////
  //
  bool Material::operator==(const Material &right) const
    {return ( (d_aeff == right.d_aeff) && 
	      (d_zeff == right.d_zeff) &&
	      (d_rho  == right.d_rho)  &&
	      (d_x0   == right.d_x0) );
    }
  
  ////////////////////////////////////////////////////////////////////////
  //
  bool Material::operator!=(const Material &right) const
    {return ( !(*this == right) );
    }
  
  
  ////////////////////////////////////////////////////////////////////////
  //
  ostream& operator<<(ostream& os, const Material* material)
    { 
      os.setf(ios::fixed,ios::floatfield);   //  Control format
      int oldprec = os.precision();          //  Save current precision

      os
	<< "Material: " << material->d_name 
	<< " (" << material->d_symbol << ")"
	<< ", Z = " << setprecision(2) << material->d_zeff 
	<< ", A = " << setprecision(2) << material->d_aeff
	<< ", density = " << setprecision(2) << material->d_rho << " g/cm**3"
	<< ", X0 = " << setprecision(2) << material->d_x0 << " g/cm**2";
      
      os.unsetf(ios::fixed);     //  Reset format       
      os.precision(oldprec);     //  Reset precision

      return os;
    }
  
  
  ////////////////////////////////////////////////////////////////////////
  //
  ostream& operator<<(ostream& os, const Material& material)
    {
      os << &material;        
      return os;
    }
  
}  //  namespace D0Material