#!/bin/python
from numpy import loadtxt, interp, zeros, size, savetxt, log, trapz, array, pi, nan_to_num 
from src.read import select_events, ReadProfile
from src.distribution import distribution
from src.analytic import distance, Mpc, MeV

EBL_models=["Dominguez"]
redshifts=["0.001168203","0.001561193","0.002086388","0.002788261","0.003726248", 
   "0.004979780","0.006655007","0.008893791","0.011885710","0.015884140", 
   "0.021227660","0.028368760","0.037912180","0.050666060","0.067710430",
   "0.090488610","0.120929500","0.161610900","0.215977800","0.288634000",
   "0.385732100","0.515494700","0.688910100","0.920663500","1.230380000",
   "1.644287000","2.197435000","2.936665000","3.924576000","5.244827000"]  

Emin = 0.1074077 # GeV
nE= 50

SaveResults=zeros((size(redshifts)+1,nE+1))    
SaveResults[0,0] = 0
SaveResults0=zeros((size(redshifts)+1,nE+1))    
SaveResults0[0,0] = 0

# Ajello et Al matrix
# > Emissivity (ph/MeV/cm^2/s/sr/bin-z) as a function of redshift and Observed Energy (MeV)
# > First line are Energies in MeV
# > First element of each subsequent line is the REDSHIFT
mat = loadtxt("Results/EmissMatrix__PLE_Paper_10TeV.txt")
Eajello = mat[0,1:]*1e-3 # MeV --> GeV
Emiss = nan_to_num(mat[1:,1:])*1e3   # ph/MeV/cm^2/s/sr/bin-z --> ph/GeV/cm^2/s/sr/bin-z 

for EBL in EBL_models:
   print "EBL model:", EBL 
   i=1
   for z in redshifts:
      fileId="Simulations/z="+z+"/"+EBL
      Emax = ReadProfile(fileId,cols=[1])*1e3 # GeV
      print "\t > redshift", z, "Emax =", Emax, "GeV", log(Emax/Emin)
      weight, energy, time, theta, phi, Esource, gen = select_events(fileId,Erange=[Emin,Emax])
      weight *= Esource*log(Emax/Emin) # such as int_Emin^Emax dNs/dE dE = 1
      # recompute weight according to Ajello et Al matrix
      dNdE_ajello = interp(Esource,Eajello*(1+float(z)),Emiss[i-1,:])/(1+float(z))
      weight *= dNdE_ajello#*16*pi**2*(distance(float(z))[1]*Mpc)**2
      E,dNdE = distribution(energy,weight,nbBins=nE,x_range=[Emin,110e3])   
      SaveResults[0,1:] = E
      SaveResults[i,1:] = dNdE
      SaveResults[i,0] = float(z)
      cond = (gen==0)
      E,dNdE = distribution(energy[cond],weight[cond],nbBins=nE,x_range=[Emin,110e3])   
      SaveResults0[0,1:] = E
      SaveResults0[i,1:] = dNdE
      SaveResults0[i,0] = float(z)
      
      i+=1

   saving_file = "Results/EGB_"+EBL+"-gen=0.dat"   
   savetxt(saving_file,SaveResults0,fmt='%1.4e')
   saving_file = "Results/EGB_"+EBL+".dat"   
   savetxt(saving_file,SaveResults,fmt='%1.4e')
   print "\t >>> Results saved in ", saving_file