from numpy import log10, histogram
from matplotlib.pyplot import figure, show
from matplotlib import gridspec
from Read import ReadRadial, ReadAngleVsEnergy
from Analytic import Analytic_theta
from Constants import degre

Eliminf = 1e0 #GeV
Elimsup = 1e5 #GeV

def angle_vs_energy(theta,energy,weight,nbPhotonsEmitted=1,nbBins=40):
   '''
      Compute arrival angle versus energy
      Input: directory name
      Output: energy, arrival angle
   '''
   energy =log10(energy)
   angle,ener =histogram(energy,nbBins,weights=weight*theta)
   nb,ener =histogram(energy,nbBins,weights=weight)
   ener = 10**ener
   enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2
   angle /= nb * nbPhotonsEmitted

   return enercenter, angle

def drawAngle_vs_energy(files):
   '''
      Plot arrival angle versus energy + analytic expression
      Input: list of directories
      Output: graph of arrival angle versus energy
   '''
   fig = figure()
   ax1 = fig.add_subplot(111)

   for fileId in files:
      ener,angle = ReadAngleVsEnergy(fileId,[0,1])
      p=ax1.plot(ener,angle,drawstyle='steps-mid',label=fileId)

      yfit = Analytic_theta(ener,fileId)
      ax1.plot(ener,yfit,'--',color="m",linewidth=2,label="analytic")

   ax1.legend(loc="best")
   ax1.set_xscale('log')
   ax1.set_yscale('log')
   ax1.set_ylim([0,180])
   ax1.grid(b=True,which='major')
   ax1.set_ylabel("$\\theta$ [rad]")
   ax1.set_xlabel("energy [GeV]")

   show()