from numpy import histogram, log10, shape
from matplotlib.pyplot import figure, show
from scipy.integrate import quad
from Read import ReadTime, ReadWeight, Readz_source, ReadNphot_source
from Constants import yr
from Integrand import comobileTime

def timing(fileName):
   time=ReadTime(fileName)
   weight=ReadWeight(fileName)
   prop = float(shape(time[time<0])[0])/shape(time)[0]
   print "file", fileName, "->", shape(time)[0], "events", "negative time:", shape(time[time<0])[0], "~", prop

   time=time[time>0]
   weight=weight[time>0]

   nbBins = 100
   zSource = Readz_source(fileName)
   tlim = quad(comobileTime,zSource,0)[0]/yr
   time = log10(time/tlim)

   dN,dt=histogram(time,nbBins,weights=weight)
   timecenter=(dt[1:nbBins+1]+dt[0:nbBins])/2
   binSize=dt[1:nbBins+1]-dt[0:nbBins]
   Nmax=ReadNphot_source(fileName) # Nb photons emitted by the source
   dNdt=dN/(Nmax*binSize)
   maxflux = ReadNphot_source(fileName)
   dNdt /= maxflux

   return timecenter, dNdt

def drawTiming(files):
   fig = figure()
   ax = fig.add_subplot(111)
   for fileId in files:
      time,dNdt = timing(fileId) 
      ax.plot(time[dNdt!=0],dNdt[dNdt!=0],linestyle="steps-mid",label=fileId)

   ax.set_yscale('log')
   ax.grid(b=True,which='major')
   ax.legend(loc="best")
   ax.set_xlabel("Time [$\Delta t / t$ log scale]")
   ax.set_ylabel("$t\ dN/dt$")

   show()