// requires ROOT: http://root.cern.ch/
//
//    usage:            compile using ROOT by: ".L library_Xo.C+" (note the plus sign)
//
//      By:             Jaideep Singh
//
//      Date:           03/07/07
//
//      descp:          Radiation Thickness Calculator
//
//                      To use this library, place the following line in your macro to load the lib:
//
//                              gSystem->Load("library_Xo_C.so");
//
//                      Make sure that it has been compiled as described in the usage above.
//                      Make sure that the path is correct in the "Load" command
//                      Make sure that you also have the "library_atomic.C" as well
//
//                      calculated radiation thickness of elements and the following composite materials:
//                          H2, D2, CO2, SiO2, Na2O, CaO, BaO, SrO, BaO, B2O, K2O, Al2O3, As2O3, 
//                          air, kapton, polystyrene, C1720, GE180, H20
//
//      references:     Davies, et al, ,Phys. Rev. 93, 788-95 (1954)
//                      Tsai, Rev. Mod. Phys. 46, 815-51 (1974), his Eq. 3.3 for f(Z) is wrong
//                      Erratum: Tsai Rev. Mod. Phys. 49, 421 (1977)
//                      Particle Data Group, "Review of Particle Physics" Phys. Lett. B Vol 592 July 2004
//                      Material Composition Data (NIST 03/07/07): http://physics.nist.gov/cgi-bin/Star/compos.pl


// head files needed to compile macro
#include "TMath.h"
#include "Riostream.h"
#include "library_pc.C"
#include "library_atomic.C"

// Coulomb Correction TSAI74(3.3)
// from: "Theory of Bremsstrahlung and Pair Production. II. Integral Cross Section for Pair Production"
//   by: Handel Davies, H.A. Bethe, and L.C. Maximon
//  ref: Phys. Rev. 93, 788-95 (1954)
// NOTE: Tsai Rev. Mod. Phys. 46, 815-51 (1974), his Eq. 3.3 for f(Z) is wrong
//       it is corrected in the subsequent Erratum: Tsai Rev. Mod. Phys. 49, 421 (1977)
Double_t fDBM(Int_t Z)
{
	Int_t nu_max = 10000;
	Double_t z = Z*pc_alpha();
	Double_t output = 0.0;
	if (Z >= 1 && Z <= 92)
	{
		for (Int_t nu = 1 ; nu < nu_max ; nu++)
		{	output += 1.0/nu/(z*z+nu*nu);	}
		output *= z*z;
	}
	else
	{ 
		output = 0.0;
		cerr << "Coulomb correction f(Z) is valid only for Z >= 1 and Z <= 92, but not Z = " << Z << endl;
	}
	return output;
}
// Radiation Logarithm TSAI74(3.67 and Table B.2)
Double_t Lrad(Int_t Z)
{
	Double_t output = TMath::Log(184.15) - 1.0/3.0*TMath::Log(Z);
	if (Z == 1)
	{	output = 5.31;	}
	else if (Z == 2)
	{	output = 4.79;	}
	else if (Z == 3)
	{	output = 4.74;	}
	else if (Z == 4)
	{	output = 4.71;	}
	return output;
}

// Radiation Logarithm prime TSAI74(3.68 and Table B.2)
Double_t Lrad_prime(Int_t Z)
{
	Double_t output = TMath::Log(1194.0) - 2.0/3.0*TMath::Log(Z);
	if (Z == 1)
	{	output = 6.144;	}
	else if (Z == 2)
	{	output = 5.621;	}
	else if (Z == 3)
	{	output = 5.805;	}
	else if (Z == 4)
	{	output = 5.924;	}
	return output;
}

// total Bremsstrahlung cross section in barns
Double_t sigma_rad(Int_t Z)
{
	Double_t scale = 4.0*pc_alpha()*pc_re()*pc_re();	// m^2
	Double_t output = Z*Z*( Lrad(Z)-fDBM(Z) ) + Z*Lrad_prime(Z);
	return scale*output/1.0E-28;	// barns
}

// Bremsstrahlung spectrum "b" scale TSAI71
Double_t fb_scale(Int_t Z)
{
	Double_t scale = 4.0/3.0;
	Double_t top = Z + 1.0;
	Double_t bot = Z*( Lrad(Z)-fDBM(Z) ) + Lrad_prime(Z);
	return scale*(1.0+top/bot/12.0);
}

// radiation length in g/cm^2 of a given isotope by Z and A
Double_t Xo(Int_t Z , Double_t A)
{
	Double_t scale = A/pc_NA(); 	// g
	Double_t output = sigma_rad(Z);	// barns
	output *= 1.0E-28;		// m^2
	output *= 100.0*100.0;		// cm^2
	return scale/output;		// g/cm^2
}

// radiation length in g/cm^2 of a given isotope by name
Double_t Xo_isotope(TString isotope)
{
	Int_t    Z = TMath::Nint(atomic(isotope,"Z"));	// atomic returns all parameters as Double_t, whereas Z is usually Int_t
	Double_t A =             atomic(isotope,"A");
	return Xo(Z,A);
}

// Bremsstrahlung spectrum "b" scale of a given isotope by name TSAI71
Double_t fb_isotope(TString isotope)
{
	Int_t    Z = TMath::Nint(atomic(isotope,"Z"));	// atomic returns all parameters as Double_t, whereas Z is usually Int_t
	return fb_scale(Z);
}

// radiation length in g/cm^2 of a composite material using Bragg's Rule
// PDG2004: Particle Data Group, "Review of Particle Physics" Phys. Lett. B Vol 592 July 2004
Double_t Xo_composite(TString material)
{	
	const Int_t max_comp = 8;
	Int_t num_comp = 0;
	Double_t Xo_array[max_comp];
	Double_t pctW[max_comp];
	if (material == "H2")	    // molecular H2 corrected for binding effects, PDG2004 p.98-9
	{
		num_comp = 1;
		Xo_array[0] = 61.28;
		pctW[0] = 1.0;
	}
	else if (material == "D2")  // molecular D2 corrected for binding effects, PDG2004 p.98-9
	{
		num_comp = 1;
		Xo_array[0] = 61.28*atomic("D","A")/atomic("H","A");
		pctW[0] = 1.0;
	}
	else if (material == "CO2")
	{
		num_comp = 2;
		Xo_array[0] = Xo_isotope("C");
		pctW[0] = atomic("C","A");
		Xo_array[1] = Xo_isotope("O");
		pctW[1] = 2.0*atomic("O","A");
	}
	else if (material == "SiO2")
	{
		num_comp = 2;
		Xo_array[0] = Xo_isotope("Si");
		pctW[0] = atomic("Si","A");
		Xo_array[1] = Xo_isotope("O");
		pctW[1] = 2.0*atomic("O","A");
	}
	else if (material == "Na2O")
	{
		num_comp = 2;
		Xo_array[0] = Xo_isotope("Na");
		pctW[0] = 2.0*atomic("Na","A");
		Xo_array[1] = Xo_isotope("O");
		pctW[1] = atomic("O","A");
	}
	else if (material == "CaO")
	{
		num_comp = 2;
		Xo_array[0] = Xo_isotope("Ca");
		pctW[0] = atomic("Ca","A");
		Xo_array[1] = Xo_isotope("O");
		pctW[1] = atomic("O","A");
	}
	else if (material == "BaO")
	{
		num_comp = 2;
		Xo_array[0] = Xo_isotope("Ba");
		pctW[0] = atomic("Ba","A");
		Xo_array[1] = Xo_isotope("O");
		pctW[1] = atomic("O","A");
	}
	else if (material == "SrO")
	{
		num_comp = 2;
		Xo_array[0] = Xo_isotope("Sr");
		pctW[0] = atomic("Sr","A");
		Xo_array[1] = Xo_isotope("O");
		pctW[1] = atomic("O","A");
	}
	else if (material == "MgO")
	{
		num_comp = 2;
		Xo_array[0] = Xo_isotope("Mg");
		pctW[0] = atomic("Mg","A");
		Xo_array[1] = Xo_isotope("O");
		pctW[1] = atomic("O","A");
	}
	else if (material == "B2O3")
	{
		num_comp = 2;
		Xo_array[0] = Xo_isotope("B");
		pctW[0] = 2.0*atomic("B","A");
		Xo_array[1] = Xo_isotope("O");
		pctW[1] = 3.0*atomic("O","A");
	}
	else if (material == "Al2O3")
	{
		num_comp = 2;
		Xo_array[0] = Xo_isotope("Al");
		pctW[0] = 2.0*atomic("Al","A");
		Xo_array[1] = Xo_isotope("O");
		pctW[1] = 3.0*atomic("O","A");
	}
	else if (material == "K2O")
	{
		num_comp = 2;
		Xo_array[0] = Xo_isotope("K");
		pctW[0] = 2.0*atomic("K","A");
		Xo_array[1] = Xo_isotope("O");
		pctW[1] = atomic("O","A");
	}
	else if (material == "As2O3")
	{
		num_comp = 2;
		Xo_array[0] = Xo_isotope("As");
		pctW[0] = 2.0*atomic("As","A");
		Xo_array[1] = Xo_isotope("O");
		pctW[1] = 3.0*atomic("O","A");
	}
	else if (material == "C1720")	// Corning 1729 glass composition from ROMALIS97, p. 97
	{
		num_comp = 8;
		Xo_array[0] = Xo_composite("SiO2");
		pctW[0] = 0.599;
		Xo_array[1] = Xo_composite("Na2O");
		pctW[1] = 0.01;
		Xo_array[2] = Xo_composite("CaO");
		pctW[2] = 0.074;
		Xo_array[3] = Xo_composite("MgO");
		pctW[3] = 0.088;
		Xo_array[4] = Xo_composite("B2O3");
		pctW[4] = 0.047;
		Xo_array[5] = Xo_composite("Al2O3");
		pctW[5] = 0.182;
		Xo_array[6] = Xo_composite("K2O");
		pctW[6] = 0.0;
		Xo_array[7] = Xo_composite("As2O3");
		pctW[7] = 0.0;
	}
	else if	(material == "GE180")	// from Al Tobais 02/15/07
	{
		num_comp = 5;
		Xo_array[0] = Xo_composite("SiO2");
		pctW[0] = 0.605;
		Xo_array[1] = Xo_composite("CaO");
		pctW[1] = 0.065;
		Xo_array[2] = Xo_composite("Al2O3");
		pctW[2] = 0.144;
		Xo_array[3] = Xo_composite("BaO");
		pctW[3] = 0.183;
		Xo_array[4] = Xo_composite("SrO");
		pctW[4] = 0.003;
	}
	else if (material == "polystyrene") // PDG2004 p.99
	{
		num_comp = 2;
		Xo_array[0] = Xo_isotope("C");
		pctW[0] = 8.0*atomic("C","A");
		Xo_array[1] = Xo_isotope("H");
		pctW[1] = 8.0*atomic("H","A");
	}
	else if (material == "kapton" || material == "polyimide film") // PDG2004 p.99
	{
		num_comp = 4;
		Xo_array[0] = Xo_isotope("C");
		pctW[0] = 22.0*atomic("C","A");
		Xo_array[1] = Xo_isotope("H");
		pctW[1] = 10.0*atomic("H","A");
		Xo_array[2] = Xo_isotope("N");
		pctW[2] = 2.0*atomic("N","A");
		Xo_array[3] = Xo_isotope("O");
		pctW[3] = 5.0*atomic("O","A");
	}
	else if (material == "H2O" || material == "ice" || material == "water")
	{
		num_comp = 2;
		Xo_array[0] = Xo_isotope("H");
		pctW[0] = 2.0*atomic("H","A");
		Xo_array[1] = Xo_isotope("O");
		pctW[1] = atomic("O","A");
	}
	else if (material == "air")  // dry air, NIST estar materials page
	{
		num_comp = 4;
		Xo_array[0] = Xo_isotope("C");
		pctW[0] = 0.000124;
		Xo_array[1] = Xo_isotope("N");
		pctW[1] = 0.755267;
		Xo_array[2] = Xo_isotope("O");
		pctW[2] = 0.231781;
		Xo_array[3] = Xo_isotope("Ar");
		pctW[3] = 0.012827;
	}
	else
	{ 	cerr << "(Xo_composite) Unknown material: " << material << endl;	}
	Double_t total = 0.0;
	Double_t output = 0.0;
	for (Int_t n = 0 ; n < num_comp ; n++)
	{
		total  += pctW[n];
		output += pctW[n]/Xo_array[n];
	}
	return total/output;
}