#!/usr/bin/python

from sys import argv
import matplotlib.pyplot as plt
from matplotlib import gridspec
from Constantes import *

if len(argv)<2:
   print "Please give at least 1 file id \n For example: 15 for Results_EGMF15"
   exit()

nbBins = 100

Eelec,Eperp,Egamma=loadtxt("Results_EGMF"+argv[1]+"/Compton_Energy",unpack=True)
deltaTheta,theta,dIC,RLtrue,RLestimate=loadtxt("Results_EGMF"+argv[1]+"/Compton_Lambda",unpack=True)

B=10**(-float(argv[1]))
lambdaB=1 #Mpc

# Distance between 2 IC versus Lepton energy 
#============================================
fig = plt.figure()
gs = gridspec.GridSpec(3, 1, height_ratios=[4,1,1]) 
fig.subplots_adjust(hspace=0.001)
ax = fig.add_subplot(gs[0])
Edata =log10(Eelec)
dist,ener =histogram(Edata,nbBins,weights=dIC)
nb,ener =histogram(Edata,nbBins)
ener = 10**ener
enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2
dist = dist/nb
ax.plot(enercenter,dist,".b",label="$\\lambda_{IC}$ - Simulation")

def fit_lambdaIC(E):
   return lambdaIC*E**0# Mpc

yfit = fit_lambdaIC(enercenter)
ax.plot(enercenter,yfit,'--c',linewidth=2,
      label="$\\lambda_{IC}= %.2f$ kpc"%(lambdaIC*1e3))

ax2 = plt.subplot(gs[1], sharex=ax)
error=dist/yfit-1
ax2.plot(enercenter,error,'+b')
step=0.2*(max(error)-min(error))
plt.ylim(min(error)-step, max(error)+step)

def fit_Dic(E):
   return Dic*(E/1000)**(-1) #Mpc

yfit = fit_Dic(enercenter)
ax.plot(enercenter,yfit,'--g',linewidth=2,
      label="$D_{IC}=%.2f \\left(\\frac{E}{1TeV}\\right)^{-1}$Mpc"%Dic)

Edata =log10(Eperp)
dist,ener =histogram(Edata,nbBins,weights=RLestimate)
nb,ener =histogram(Edata,nbBins)
ener = 10**ener
enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2
dist = dist/nb
ax.plot(enercenter,dist,".r",label="$R_{L}$ - estimate")
#ax.plot(Eperp,RLestimate,".r",label="$R_{L}$ - estimate")
#ax.plot(Eperp,RLtrue,".y",label="$R_{L}$ - true")

def fit_RL(E,B):
   return RL*(E/1e3)*(B/1e-17)**-1 # Mpc

yfit = fit_RL(enercenter,B)
ax.plot(enercenter,yfit,'--m',linewidth=2,
      label="$R_{L}=%.2f \\left(\\frac{E}{1TeV}\\right)$ Mpc"%(RL*1e-17/B))

ax3 = plt.subplot(gs[2], sharex=ax)
error=dist/yfit-1
ax3.plot(enercenter,error,'+r')
step=0.2*(max(error)-min(error))
plt.ylim(min(error)-step, max(error)+step)

ax.legend(loc="best",title="$10^{-%.0f"%float(argv[1])+"}$Gauss")
ax.grid(b=True,which='major')
ax2.grid(b=True,which='major')
ax3.grid(b=True,which='major')
ax.set_xscale('log')
ax.set_yscale('log')
ax2.set_xscale('log')
ax3.set_xscale('log')
ax3.set_xlabel("energy [GeV]")
ax.set_ylabel("$\lambda$ [Mpc]")
ax2.set_ylabel("relative error $d_{IC}$")
ax3.set_ylabel("relative error $R_L$")
xticklabels = ax.get_xticklabels()+ax2.get_xticklabels()
plt.setp(xticklabels, visible=False)

# Deflexion angle versus Lepton energy 
#======================================
fig = plt.figure()
gs = gridspec.GridSpec(2, 1, height_ratios=[4,1]) 
fig.subplots_adjust(hspace=0.001)
ax = fig.add_subplot(gs[0])

angle,ener =histogram(Edata,nbBins,weights=theta)
nb,ener =histogram(Edata,nbBins)
ener = 10**ener
enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2
angle = angle/nb*degre
ax.plot(enercenter,angle,".b",linewidth=2,label="Simulation -$10^{-%.0f"%float(argv[1])+"}$Gauss")

def fit_deltaIC(E):
   return delta*(E/1e3)**(-1)*(B/1e-17)*(lambdaB/0.3)**(1/2) #degre

yfit = fit_deltaIC(enercenter)
ax.plot(enercenter,yfit,'--c',linewidth=2,label="Analytic")

ax2 = plt.subplot(gs[1], sharex=ax)
error=angle/yfit-1
ax2.plot(enercenter,error,'+b')
step=0.2*(max(error)-min(error))
plt.ylim(min(error)-step, max(error)+step)

#a=0.2 # degre
#def fitLeptonDN13(E,a,B):
#   return a* (E/1000)**(-2) *(B/1e-17) 
#
#yfit = fitLeptonDN13(enercenter,a,B)
#ax.plot(enercenter,yfit,'--g',linewidth=2,label="Durrer 2013")

ax.legend(loc="best")
ax.grid(b=True,which='major')
ax2.grid(b=True,which='major')
ax.set_xscale('log')
ax2.set_xscale('log')
ax.set_yscale('log')
ax2.set_xlabel("energy [GeV]")
ax.set_ylabel("$\\delta$ [degre]")
ax2.set_ylabel("relative error $\\theta$")
xticklabels = ax.get_xticklabels()
plt.setp(xticklabels, visible=False)

# Histo leptons energy
#======================
#fig = plt.figure()
#ax = fig.add_subplot(111)
#
#plt.hist(Eelec, nbBins, log=True, facecolor="green", alpha=0.5, 
#      label="$E_{lepton}$")
#
#ax.set_title("mean E lepton: %.3f"%(mean(Eelec))+"GeV")
#ax.legend(loc="best")
#ax.grid(b=True,which='major')
#ax.set_ylabel("Nb events")
#ax.set_xlabel("$E$ [GeV]")

plt.show()