Rewrite Analysis.py :

> only one power law spectrum index > compute results gen by gen and by energy range when it makes sense > add observables.py to plot 3D space of arrival parameters (E, dt, theta)

Rewrite Analysis.py :
> only one power law spectrum index > compute results gen by gen and by energy range when it makes sense > add observables.py to plot 3D space of arrival parameters (E, dt, theta)
Thomas Fitoussi
1 parent a44230b7
Showing 23 changed files with 327 additions and 346 deletions Show diff stats
.gitignore
Analysis.py
EBL-spectrums.py
Modules/Analytic.pyc
Modules/Angle.py
Modules/Angle.pyc
Modules/Constants.pyc
Modules/Generation.pyc
Modules/Integrand.pyc
Modules/Map.py
Modules/Map.pyc
Modules/Observables.py
Modules/Read.py
Modules/Read.pyc
Modules/Spectrum.py
Modules/Spectrum.pyc
Modules/Timing.py
Modules/Timing.pyc
Modules/Weight.pyc
Modules/__init__.py
@@ -0,0 +1,2 @@
+*.pyc
+*.swp
@@ -10,21 +10,18 @@ from Modules.Angle import angle_vs_energy, radial
 from Modules.Timing import timing
 from Modules.Constants import degre
-def InProgress(i,Nmax):
-   print "  ", (i*100)/Nmax, "% done"
-
 if shape(argv)[0] < 2: 
    print "not enough arguments (at least 1)"
    exit()
-PowerSpectrum=[1,1.5,2,2.5] 
+powerlaw_index=2 
 for fileId in argv[1:]:
-   print "#=============================================================================#"
-   print "#    Analysis of", fileId
-   print "#=============================================================================#"
+   print "#===============================================#"
+   print "# Analysis of", fileId
+   print "#===============================================#"
    # read files
-   print "# 1. Reading data" 
+   print "   > Reading data" 
    time = ReadTime(fileId)
    energy = ReadEnergy(fileId)
    weightini, generation, theta_arrival, Esource, dir_source = ReadExtraFile(fileId,[2,3,4,5,6])
@@ -39,228 +36,135 @@ for fileId in argv[1:]:
    theta = thetaDir - thetaPos
    phi = phiDir -phiPos
-   Gen_contrib = []
-   Gen_cont = zeros((int(max(generation))+1))
-   Source = []
-   Spectrum = []
-   Timing = []
-   Angle_Energy = []
-   Radial = []
-
-   print "# 2. Computing powerlaw spectrum" 
-   #=============================================================================#
-   #  VERSUS SOURCE SPECTRUM
-   #=============================================================================#
-   for powerlaw_index in PowerSpectrum:
-      # apply source spectrum
-      weight_source = (Esource/min(Esource))**(1-powerlaw_index)
-      weight = weightini* weight_source
-
-      # GENERATION CONTRIBUTION ====================================================#
-      NbTotEvents = sum(weight)
-      if Gen_contrib==[]:
-         Gen_contrib = arange(0,max(generation)+1,1)[:,newaxis]
-
-      for gen in list(set(generation)):
-         Gen_cont[int(gen)] = sum(weight[generation==gen])/NbTotEvents *100
-      Gen_contrib = append(Gen_contrib,Gen_cont[:,newaxis],axis=1)
-
-      #  SPECTRUM (SOURCE AND MEASURED) ============================================#
-      nbBins = 100
-      # draw source spectrum
-      Es=array(list(set(Esource)))
-      Ws= (Es/min(Es))**(1-powerlaw_index)  / ratio
-      Es,Fs = spectrum(Es,Ws,nbBins=nbBins)
-      Es=Es[:,newaxis]
-      Fs=Fs[:,newaxis]
-      if Source==[]:
-         Source = Es
-         Source = append(Source,Fs,axis=1)
-      else: 
-         Source = append(Source,Fs,axis=1)
-
-      # primary gamma-rays contribution and full spectrum
-      ener,flux,flux_0 = spectrum(energy,weight,generation,nbBins)
-      ener=ener[:,newaxis]
-      flux=flux[:,newaxis]
-      flux_0=flux_0[:,newaxis]
-      if Spectrum==[]:
-         Spectrum = ener
-         Spectrum = append(Spectrum,flux,axis=1)
-         Spectrum = append(Spectrum,flux_0,axis=1)
-      else:
-         Spectrum = append(Spectrum,flux,axis=1)
-         Spectrum = append(Spectrum,flux_0,axis=1)
-
-      #  IMAGING, RADIAL DISTRIBUTION AND ANGLE VERSUS ENERGY ======================#
-      nbBins = 50
-      cond = energy>1e-3
-      ener,angle = angle_vs_energy(theta_arrival[cond],energy[cond],weight[cond],nbBins)
-      theta2,dndtheta2 = radial(theta[cond],weight[cond],nbBins)
-      theta2=theta2[:,newaxis]
-      dndtheta2=dndtheta2[:,newaxis]
-      ener=ener[:,newaxis]
-      angle=angle[:,newaxis]
-
-      if Radial==[]:
-         Radial = theta2
-         Radial = append(Radial,dndtheta2,axis=1)
-         Angle_Energy = ener
-         Angle_Energy = append(Angle_Energy,angle,axis=1)
-      else:
-         Radial = append(Radial,dndtheta2,axis=1)
-         Angle_Energy = append(Angle_Energy,angle,axis=1)
-
-
-      #  TIME DISTRIBUTION AND TIME DELAY VERSUS ANGLE =============================#
-      nbBins = 100
-      delta_t,dNdt = timing(time,weight,nbBins)
-      delta_t=delta_t[:,newaxis]
-      dNdt=dNdt[:,newaxis]
-      if Timing==[]:
-         Timing = delta_t
-         Timing = append(Timing,dNdt,axis=1)
-      else:
-         Timing = append(Timing,dNdt,axis=1)
-
-      InProgress(PowerSpectrum.index(powerlaw_index)+1,shape(PowerSpectrum)[0])
+   powerlaw_index = 2
    #=============================================================================#
-   #  VERSUS JET OPENING ANGLE
+   #  NO SELECTION
    #=============================================================================#
-   Jet_Opening=[180,60,30] # degre (180 <=> isotrop)
-   powerlaw_index = 2
-   Elim = 1e-1 # GeV
-   thetalim = 20 # degre
-
-   # apply source spectrum
+   print "   > No selection of events ..." 
+   #  SPECTRUM (SOURCE) =========================================================#
    weight_source = (Esource/min(Esource))**(1-powerlaw_index)
    weight = weightini* weight_source
-
-   print "# 3. compute images and radial distribution" 
-   for jet_opening_angle in Jet_Opening:
-      # apply selection ( /!\ USE DECREASING VALUES OF jet_opening_angle )
-      cond = (dir_source*degre <= jet_opening_angle) #& (energy>Elim) #& (abs(theta)<thetalim) & (abs(phi)<thetalim)
-      dir_source = dir_source[cond]
-      theta_arrival = theta_arrival[cond]
-      theta = theta[cond]
-      phi = phi[cond]
-      weight = weight[cond]
-      energy = energy[cond]
-      time = time[cond]
-      generation = generation[cond]
-
-      # GENERATION CONTRIBUTION ====================================================#
-      NbTotEvents = sum(weight)
-      if Gen_contrib==[]:
-         Gen_contrib = arange(0,max(generation)+1,1)[:,newaxis]
-
-      for gen in list(set(generation)):
-         Gen_cont[int(gen)] = sum(weight[generation==gen])/NbTotEvents *100
-      Gen_contrib = append(Gen_contrib,Gen_cont[:,newaxis],axis=1)
-
-      #  SPECTRUM (SOURCE AND MEASURED) ============================================#
-      nbBins = 100
-      # draw source spectrum
-      Es=array(list(set(Esource)))
-      Ws= (Es/min(Es))**(1-powerlaw_index) / ratio
-      Es,Fs = spectrum(Es,Ws,nbBins=nbBins)
-      Es=Es[:,newaxis]
-      Fs=Fs[:,newaxis]
-      if Source==[]:
-         Source = Es
-         Source = append(Source,Fs,axis=1)
-      else: 
-         Source = append(Source,Fs,axis=1)
-
-      # primary gamma-rays contribution and full spectrum
-      ener,flux,flux_0 = spectrum(energy,weight,generation,nbBins)
-      ener=ener[:,newaxis]
-      flux=flux[:,newaxis]
-      flux_0=flux_0[:,newaxis]
-      if Spectrum==[]:
-         Spectrum = ener
-         Spectrum = append(Spectrum,flux,axis=1)
-         Spectrum = append(Spectrum,flux_0,axis=1)
-      else:
-         Spectrum = append(Spectrum,flux,axis=1)
-         Spectrum = append(Spectrum,flux_0,axis=1)
-
-      #  IMAGING, RADIAL DISTRIBUTION AND ANGLE VERSUS ENERGY ======================#
-      nbBins = 50
-      computeMap(theta,phi,weight,energy,fileId,jet_opening_angle,nbBins,borne=[thetalim,thetalim])
-      ener,angle = angle_vs_energy(theta_arrival,energy,weight,nbBins)
-      theta2,dndtheta2 = radial(theta_arrival,weight,nbBins)
-      theta2=theta2[:,newaxis]
-      dndtheta2=dndtheta2[:,newaxis]
-      ener=ener[:,newaxis]
-      angle=angle[:,newaxis]
-
-      if Radial==[]:
-         Radial = theta2
-         Radial = append(Radial,dndtheta2,axis=1)
-         Angle_Energy = ener
-         Angle_Energy = append(Angle_Energy,angle,axis=1)
-      else:
-         Radial = append(Radial,dndtheta2,axis=1)
-         Angle_Energy = append(Angle_Energy,angle,axis=1)
-
-      #  TIME DISTRIBUTION AND TIME DELAY VERSUS ANGLE =============================#
-      nbBins = 100
-      delta_t,dNdt = timing(time,weight,nbBins)
-      delta_t=delta_t[:,newaxis]
-      dNdt=dNdt[:,newaxis]
-      Timing = append(Timing,dNdt,axis=1)
-
-      InProgress(Jet_Opening.index(jet_opening_angle)+1,shape(Jet_Opening)[0])
+   nbBins = 100
+   Es=array(list(set(Esource)))
+   Ws= (Es/min(Es))**(1-powerlaw_index)  / ratio
+   Es,Fs = spectrum(Es,Ws,nbBins=nbBins)
+   Es=Es[:,newaxis]
+   Fs=Fs[:,newaxis]
+   Source = Es
+   Source = append(Source,Fs,axis=1)
+
+   #  IMAGING ===================================================================#
+   print "     ... Computing image"
+   nbBins = 50
+   thetalim = 20
+   computeMap(theta,phi,weight,energy,fileId,nbBins,borne=[thetalim,thetalim])
+
+   #  SPECTRUM (MEASURED) =======================================================#
+   print "     ... Computing spectrum"
+   nbBins = 100
+   ener,flux = spectrum(energy,weight,nbBins)
+   ener=ener[:,newaxis]
+   flux=flux[:,newaxis]
+   Spectrum = ener
+   Spectrum = append(Spectrum,flux,axis=1)
+
+   #  ANGLE VERSUS ENERGY =======================================================#
+   print "     ... Computing arrival angle versus energy"
+   nbBins = 100
+   ener,angle = angle_vs_energy(theta_arrival,energy,weight,nbBins)
+   ener=ener[:,newaxis]
+   angle=angle[:,newaxis]
+   Angle_Energy = ener
+   Angle_Energy = append(Angle_Energy,angle,axis=1)
+
+   #  RADIAL ====================================================================#
+   print "     ... Computing radial distribution"
+   nbBins = 100
+   theta2,dndtheta = radial(theta,weight,nbBins)
+   theta2=theta2[:,newaxis]
+   dndtheta=dndtheta[:,newaxis]
+   Radial = theta2
+   Radial = append(Radial,dndtheta,axis=1)
+
+   #  TIME DISTRIBUTION AND TIME DELAY VERSUS ANGLE =============================#
+   print "     ... Computing time distribution"
+   nbBins = 100
+   delta_t,dNdt = timing(time,weight,nbBins)
+   delta_t=delta_t[:,newaxis]
+   dNdt=dNdt[:,newaxis]
+   Timing = delta_t
+   Timing = append(Timing,dNdt,axis=1)
    #=============================================================================#
    #  BY ENERGY RANGE
    #=============================================================================#
-   print "# 3. By energy range" 
-   powerlaw_index = 2
-   weight_source = (Esource/min(Esource))**(1-powerlaw_index)
-   weight = weightini* weight_source
-
+   print "   > Selection by energy range ..." 
    Emin = [1e-3,1e0,1e3] #GeV
    Emax = [1e0,1e3,1e5] #GeV
+
    for n in arange(0,3,1):
       cond= (energy>Emin[n]) & (energy<Emax[n])
-      #  TIME DISTRIBUTION AND TIME DELAY VERSUS ANGLE =============================#
+      #  RADIAL =================================================================#
+      nbBins = 100
+      theta2,dndtheta = radial(theta[cond],weight[cond],nbBins)
+      dndtheta=dndtheta[:,newaxis]
+      Radial = append(Radial,dndtheta,axis=1)
+
+      #  TIME DISTRIBUTION AND TIME DELAY VERSUS ANGLE ==========================#
       nbBins = 100
       delta_t,dNdt = timing(time[cond],weight[cond],nbBins)
-      delta_t=delta_t[:,newaxis]
       dNdt=dNdt[:,newaxis]
       Timing = append(Timing,dNdt,axis=1)
+      print "     ... ", ((n+1)*100)/3,"% done"
+
    #=============================================================================#
    #  BY GENERATION
    #=============================================================================#
-   print "# 3. By generation" 
+   print "   > Selection by generation ..." 
+   NbTotEvents = sum(weight)
    gen_tab =list(set(generation))
-   powerlaw_index = 2
-   weight_source = (Esource/min(Esource))**(1-powerlaw_index)
-   weight = weightini* weight_source
-
+   Gen_cont = zeros([shape(gen_tab)[0],2])
    for gen in gen_tab:
       cond = generation==gen
-      print shape(time[cond])
+      contrib =sum(weight[cond])/NbTotEvents*100 
+      i = gen_tab.index(gen)
+      Gen_cont[i,0]=int(gen)
+      Gen_cont[i,1]=contrib
+      print "     ... gen=",int(gen),"-> contribution:",int(contrib),"%"
+
+      #  SPECTRUM (MEASURED) ====================================================#
+      nbBins = 100
+      ener,flux = spectrum(energy[cond],weight[cond],nbBins)
+      flux=flux[:,newaxis]
+      Spectrum = append(Spectrum,flux,axis=1)
-      #  TIME DISTRIBUTION AND TIME DELAY VERSUS ANGLE =============================#
+      #  ANGLE VERSUS ENERGY ====================================================#
+      nbBins = 100
+      ener,angle = angle_vs_energy(theta_arrival[cond],energy[cond],weight[cond],nbBins)
+      angle=angle[:,newaxis]
+      Angle_Energy = append(Angle_Energy,angle,axis=1)
+
+      #  RADIAL =================================================================#
+      nbBins = 100
+      theta2,dndtheta = radial(theta[cond],weight[cond],nbBins)
+      dndtheta=dndtheta[:,newaxis]
+      Radial = append(Radial,dndtheta,axis=1)
+
+      #  TIME DISTRIBUTION AND TIME DELAY VERSUS ANGLE ==========================#
       nbBins = 100
       delta_t,dNdt = timing(time[cond],weight[cond],nbBins)
-      delta_t=delta_t[:,newaxis]
       dNdt=dNdt[:,newaxis]
       Timing = append(Timing,dNdt,axis=1)
    #=============================================================================#
-   print "# 4. writing files" 
-   savetxt(resultDirectory+fileId+"/Generation.txt",Gen_contrib)
+   print "   > writing files" 
    savetxt(resultDirectory+fileId+"/Spectrum.txt",Spectrum)
    savetxt(resultDirectory+fileId+"/Source_spectrum.txt",Source)
    savetxt(resultDirectory+fileId+"/Angle_vs_Energy.txt",Angle_Energy)
    savetxt(resultDirectory+fileId+"/Radial_distribution.txt",Radial)
    savetxt(resultDirectory+fileId+"/Timing.txt",Timing)
-   print "#=============================================================================#"
-
+   savetxt(resultDirectory+fileId+"/Generation.txt",Gen_cont)
+   print "#===============================================#"
@@ -58,7 +58,7 @@ ax.plot(hv,nCMB(hv,z)*hv**2,&quot;-k&quot;,label=&quot;CMB&quot;)
 ax.set_xscale('log')
 ax.set_yscale('log')
 ax.grid(b=True,which='major')
-ax.legend(loc="best",frameon=False,framealpha=0.5)
+ax.legend(loc="best")#,frameon=False,framealpha=0.5)
 ax.set_xlabel("energy [eV]")
 ax.set_ylabel("$n$ [photon.eV.cm$^{-3}$]")
@@ -80,7 +80,7 @@ for z_index in [1,9,12,15,17]:
 ax2.set_xscale('log')
 ax2.set_yscale('log')
 ax2.grid(b=True,which='major')
-ax2.legend(loc="best",frameon=False,framealpha=0.5)
+ax2.legend(loc="best")#,frameon=False,framealpha=0.5)
 ax2.set_xlabel("energy [eV]")
 ax2.set_ylabel("$n$ [photon.eV.cm$^{-3}$]")
-from numpy import log10, histogram
+from numpy import log10, histogram, arange
 from matplotlib.pyplot import figure, show
 from matplotlib import gridspec
-from Read import ReadRadial, ReadAngleVsEnergy
+from Read import ReadRadial, ReadAngleVsEnergy, ReadGeneration
 from Analytic import Analytic_theta
 from Constants import degre
@@ -21,7 +21,7 @@ def angle_vs_energy(theta,energy,weight,nbBins=50):
    return enercenter, angle
 def radial(theta,weight,nbBins=50):
-   cond = theta!=0
+   cond = theta>0
    theta = log10(theta[cond])
    weight = weight[cond]
    dn,angle2 = histogram(theta,nbBins,range=[log10(1e-3),log10(25)],weights=weight)
@@ -35,7 +35,7 @@ def radial(theta,weight,nbBins=50):
    return thetacenter, dndtheta
-def drawRadial(files,all_source_spectrum=False,all_jets=False):
+def drawRadial(files,plot_gen=False,plot_energy_range=False):
    '''
       Plot flux versus arrival angle
       Input: list of directories
@@ -45,21 +45,24 @@ def drawRadial(files,all_source_spectrum=False,all_jets=False):
    ax = fig.add_subplot(111)
    for fileId in files:
-      if all_source_spectrum:
-         i=1
-         for powerlaw_index in [1,1.5,2,2.5]: 
+      theta,dndtheta2 = ReadRadial(fileId,[0,1])
+      ax.plot(theta,dndtheta2,drawstyle='steps-mid')
+
+      if plot_energy_range:
+         labels = ["MeV band","GeV band","TeV band"]
+         i = 2
+         for n in arange(0,3,1):
             theta,dndtheta2 = ReadRadial(fileId,[0,i])
-            ax.plot(theta,dndtheta2,drawstyle='steps-mid',label="p=%.1f"%powerlaw_index)
+            ax.plot(theta,dndtheta2,drawstyle='steps-mid',label=labels[n])
             i+=1
-      elif all_jets:
+
+      if plot_gen:
+         generation = ReadGeneration(fileId,[0])
          i=5
-         for jet_opening in ["isotrop","60 deg","30 deg"]: 
+         for gen in generation:
             theta,dndtheta2 = ReadRadial(fileId,[0,i])
-            ax.plot(theta,dndtheta2,drawstyle='steps-mid',label=jet_opening)
+            ax.plot(theta,dndtheta2,drawstyle='steps-mid',label="gen=%.0f"%gen)
             i+=1
-      else:
-         theta,dndtheta2 = ReadRadial(fileId,[0,5])
-         ax.plot(theta,dndtheta2,drawstyle='steps-mid',label=fileId)
    ax.legend(loc="best")
    ax.set_xscale('log')
@@ -70,7 +73,7 @@ def drawRadial(files,all_source_spectrum=False,all_jets=False):
    show()
-def drawAngle_vs_energy(files,all_source_spectrum=False,all_jets=False,PlotAnalytic=False):
+def drawAngle_vs_energy(files,plot_gen=False,plot_analytic=False):
    '''
       Plot arrival angle versus energy + analytic expression
       Input: list of directories
@@ -80,21 +83,18 @@ def drawAngle_vs_energy(files,all_source_spectrum=False,all_jets=False,PlotAnaly
    ax1 = fig.add_subplot(111)
    for fileId in files:
-      if all_source_spectrum:
-         i=1
-         for powerlaw_index in [1,1.5,2,2.5]: 
-            ener,angle = ReadAngleVsEnergy(fileId,[0,i])
-            p=ax1.plot(ener,angle,drawstyle='steps-mid',label="p=%.1f"%powerlaw_index)
-      elif all_jets:
-         i=5
-         for jet_opening in ["isotrop","60 deg","30 deg"]: 
+      ener,angle = ReadAngleVsEnergy(fileId,[0,1])
+      p=ax1.plot(ener,angle,drawstyle='steps-mid')
+
+      if plot_gen:
+         generation = ReadGeneration(fileId,[0])
+         i=2
+         for gen in generation:
             ener,angle = ReadAngleVsEnergy(fileId,[0,i])
-            p=ax1.plot(ener,angle,drawstyle='steps-mid',label=jet_opening)
-      else:
-         ener,angle = ReadAngleVsEnergy(fileId,[0,1])
-         p=ax1.plot(ener,angle,drawstyle='steps-mid',label=fileId)
+            ax1.plot(ener,angle,drawstyle='steps-mid',label="gen=%.0f"%gen)
+            i+=1
-      if PlotAnalytic:
+      if plot_analytic:
          yfit = Analytic_theta(ener,fileId)
          ax1.plot(ener,yfit,'--',color=p[0].get_color(),linewidth=2)
@@ -15,7 +15,7 @@ philim=float(philim)
 thetalim=float(thetalim)
 limits = [[-thetalim,thetalim],[-philim,philim]]
-def computeMap(theta,phi,weight,energy,saveId,jet_opening_angle=180,nbBins=50,borne=[180,180]):
+def computeMap(theta,phi,weight,energy,saveId,nbBins=50,borne=[180,180]):
    '''
       Plot 2D histogram of arrival direction of the photons with energy upper than Elim
       Input: list of directories
@@ -41,12 +41,8 @@ def computeMap(theta,phi,weight,energy,saveId,jet_opening_angle=180,nbBins=50,bo
    ax.set_xlabel("$\\theta$ $\\cos(\\phi)$ [deg]")
    ax.set_ylabel("$\\theta$ $\\sin(\\phi)$ [deg]")
    ax.grid(b=True,which='major')
-   if jet_opening_angle == 180:
-      ax.set_title(saveId+" - isotrop")
-      savefig(resultDirectory+saveId+"/map_isotrop.png")
-   else:
-      ax.set_title(saveId+" - jet opening angle: %d degre"%(int(jet_opening_angle)))
-      savefig(resultDirectory+saveId+"/map_"+str(int(jet_opening_angle))+".png")
+   ax.set_title(saveId+" - isotrop")
+   savefig(resultDirectory+saveId+"/map_isotrop.png")
 #def isotrop():
 #   return array([[1]]), "isotropic"
@@ -0,0 +1,56 @@
+from numpy import sort, log10
+from matplotlib.pyplot import figure, show
+from mpl_toolkits.mplot3d import Axes3D
+from Modules.Read import ReadEnergy, ReadTime, ReadExtraFile
+from Modules.Constants import degre
+
+def drawObservablesSpace(fileId):
+   fig = figure()
+   ax = fig.add_subplot(232,projection='3d')
+   ax1 = fig.add_subplot(233)
+   ax2 = fig.add_subplot(231)
+   ax3 = fig.add_subplot(235)
+
+   energy = ReadEnergy(fileId)
+   delay  = ReadTime(fileId)
+   generation, theta = ReadExtraFile(fileId,[3,4])
+   theta*=degre
+
+   colors=['b','g','r','m','c','y']
+   i=0
+
+   for gen in sort(list(set(generation))):
+      cond = (generation==gen) & (delay>0) & (theta>0)
+      ax.scatter(log10(energy[cond]),log10(theta[cond]),log10(delay[cond]),color=colors[i],label="gen = %.0f"%gen)
+      ax1.scatter(energy[cond],delay[cond],color=colors[i])
+      ax2.scatter(theta[cond],delay[cond],color=colors[i])
+      ax3.scatter(theta[cond],energy[cond],color=colors[i])
+      i+=1
+
+   ax.legend(loc="best")
+   #ax.set_xscale('log')
+   #ax.set_yscale('log')
+   #ax.set_zscale('log')
+   ax.set_xlabel("energy [GeV]")
+   ax.set_ylabel("$\\theta$ [deg]")
+   ax.set_zlabel("Time delay [s]")
+
+   ax1.grid(b=True,which='major')
+   ax1.set_xscale('log')
+   ax1.set_yscale('log')
+   ax1.set_xlabel("energy [GeV]")
+   ax1.set_ylabel("Time delay [s]")
+
+   ax2.grid(b=True,which='major')
+   ax2.set_xscale('log')
+   ax2.set_yscale('log')
+   ax2.set_xlabel("$\\theta$ [deg]")
+   ax2.set_ylabel("Time delay [s]")
+
+   ax3.grid(b=True,which='major')
+   ax3.set_xscale('log')
+   ax3.set_yscale('log')
+   ax3.set_xlabel("$\\theta$ [deg]")
+   ax3.set_ylabel("energy [GeV]")
+
+   show()
@@ -60,5 +60,5 @@ def ReadRadial(fileName,cols=[0,1,2,3,4]):
 def ReadTiming(fileName,cols=[0,1,2,3,4]): 
    return loadtxt(resultDirectory+fileName+"/Timing.txt",unpack=True,usecols=cols)
-def ReadGeneration(fileName,cols=[0,1,2,3,4]):
+def ReadGeneration(fileName,cols=[0,1]):
    return loadtxt(resultDirectory+fileName+"/Generation.txt",unpack=True,usecols=cols)
 from numpy import histogram, log10, sqrt, array
 from matplotlib.pyplot import figure, show, setp
 from matplotlib import gridspec
-from Read import ReadSpectrum, ReadSourceSpectrum
+from Read import ReadSpectrum, ReadSourceSpectrum, ReadGeneration
 from Constants import degre
-from Analytic import Ethreshold_gg
+from Analytic import Ethreshold_gg, Analytic_flux, Eic
-def spectrum(energy,weight,generation=[],nbBins=100):
+def spectrum(energy,weight,nbBins=100):
    '''
       Compute flux versus energy
       Input: events energy and weight (optional: number of bins)
@@ -19,20 +19,9 @@ def spectrum(energy,weight,generation=[],nbBins=100):
    enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2
    binSize=ener[1:nbBins+1]-ener[0:nbBins]
    flux = enercenter**2 * flux/binSize 
+   return enercenter,flux
-   if generation == []:
-      return enercenter,flux
-   else:
-      cond=(generation==0)
-      flux_0,ener_0=histogram(energy[cond],nbBins,range=[log10(Emin),log10(Emax)],weights=weight[cond])
-      ener_0 = 10**ener_0
-      enercenter_0=(ener_0[1:nbBins+1]+ener_0[0:nbBins])/2
-      binSize=ener_0[1:nbBins+1]-ener_0[0:nbBins]
-      flux_0 = enercenter_0**2 * flux_0/binSize 
-      return enercenter,flux,flux_0
-
-
-def drawSpectrum(files,all_source_spectrum=False,all_jets=False,PlotAnalytic=False):
+def drawSpectrum(files,plot_gen=False,plot_source=False,plot_analytic=False):
    '''
       Plot flux versus energy
       Input: list of directories
@@ -47,59 +36,36 @@ def drawSpectrum(files,all_source_spectrum=False,all_jets=False,PlotAnalytic=Fal
    ax1 = fig.add_subplot(111)
    for fileId in files:
-      if all_source_spectrum:
-         i=1
-         j=1
-         for powerlaw_index in [1,1.5,2,2.5]: 
-            # read files
-            Es,Fs = ReadSourceSpectrum(fileId,[0,i])
-            p = ax1.plot(Es,Fs,linestyle=':')
-            # draw full spectrum
-            ener,flux,flux_0 = ReadSpectrum(fileId,[0,j,j+1])
-            ax1.plot(ener,flux,color=p[0].get_color(),drawstyle='steps-mid',label="p=%.1f"%powerlaw_index)
-            # primary gamma-rays contribution
-            #ax1.plot(ener,flux_0,color=p[0].get_color(),linestyle='--',drawstyle='steps-mid')
-            #contrib=flux_0/flux
-            #ax2.plot(ener,contrib,drawstyle='steps-mid',color=p[0].get_color())
-            i+=1
-            j+=2
-      elif all_jets:
-         i=5
-         j=9
-         for jet_opening in ["isotrop","60 deg","30 deg"]: 
+      # draw full spectrum
+      ener,flux = ReadSpectrum(fileId,[0,1])
+      p = ax1.plot(ener,flux,drawstyle='steps-mid')
+
+      if plot_gen:
+         generation = ReadGeneration(fileId,[0])
+         j=2
+         for gen in generation:
             # read files
-            Es,Fs = ReadSourceSpectrum(fileId,[0,i])
-            p = ax1.plot(Es,Fs,linestyle=':')
-            # draw full spectrum
-            ener,flux,flux_0 = ReadSpectrum(fileId,[0,j,j+1])
-            ax1.plot(ener,flux,color=p[0].get_color(),drawstyle='steps-mid',label=jet_opening)
-            # primary gamma-rays contribution
-            #ax1.plot(ener,flux_0,color=p[0].get_color(),linestyle='--',drawstyle='steps-mid')
-            #contrib=flux_0/flux
-            #ax2.plot(ener,contrib,drawstyle='steps-mid',color=p[0].get_color())
-            i+=1
-            j+=2
+            ener,flux = ReadSpectrum(fileId,[0,j])
+            ax1.plot(ener,flux,drawstyle='steps-mid',label="gen=%.0f"%gen)
+            j+=1
-      else:
-         # read files
-         #Es,Fs = ReadSourceSpectrum(fileId,[0,3])
-         #p = ax1.plot(Es,Fs,linestyle=':')
-         # draw full spectrum
-         ener,flux,flux_0 = ReadSpectrum(fileId,[0,5,6])
-         ax1.plot(ener,flux,drawstyle='steps-mid')
-         #ax1.plot(ener,flux,color=p[0].get_color(),drawstyle='steps-mid',label=fileId)
-         # primary gamma-rays contribution
-         #ax1.plot(ener,flux_0,color=p[0].get_color(),linestyle='--',drawstyle='steps-mid')
-         #contrib=flux_0/flux
-         #ax2.plot(ener,contrib,drawstyle='steps-mid',color=p[0].get_color())
+      if plot_source:
+         Es,Fs = ReadSourceSpectrum(fileId,[0,1])
+         ax1.plot(Es,Fs,color=p[0].get_color(),linestyle=':')
-   if PlotAnalytic:
+   if plot_analytic:
       minEner=min(ener)
-      alpha=flux[ener==minEner]/sqrt(minEner)
-      ax1.plot(ener,alpha*sqrt(ener),'--m',linewidth=2)
-      ax1.plot(ener,ener**(0)*max(flux)*0.95,'--m',linewidth=2)
-      ax1.axvline(x=Ethreshold_gg(), ymin=0., ymax = 1., color='r', linewidth=2)
-      ax1.axvline(x=1e4, ymin=0., ymax = 1., color='g', linewidth=2)
+      #alpha=flux[ener==minEner]/sqrt(minEner)
+      #ax1.plot(ener,alpha*sqrt(ener/100),'.-g',linewidth=2)
+      ax1.plot(ener,2*Analytic_flux(ener),'--g',linewidth=1)
+      Elim= Eic(Ethreshold_gg()/2)
+      alpha = 2*4.63e3/Elim**(1/4.)
+      ax1.plot(ener,alpha*sqrt(ener),'--r',linewidth=1)
+      alpha = 2*14.6e3
+      ax1.plot(ener,alpha*(ener/100)**(1/4.),'--m',linewidth=1)
+      #ax1.plot(ener,ener**(0)*max(flux)*0.95,'--m',linewidth=2)
+      #ax1.axvline(x=Ethreshold_gg(), ymin=0., ymax = 1., color='r', linewidth=2)
+      #ax1.axvline(x=1e4, ymin=0., ymax = 1., color='g', linewidth=2)
    ax1.legend(loc="best")
    ax1.set_xscale('log')
-from numpy import histogram, log10, shape, arange, pi, convolve
+from numpy import histogram, log10, shape, arange, pi, convolve, r_
 from matplotlib.pyplot import figure, show
-from Read import ReadTime, ReadTiming, ReadExtraFile, ReadProfile, ReadEnergy
+from Read import ReadTime, ReadTiming, ReadExtraFile, ReadProfile, ReadEnergy, ReadGeneration
 from Constants import yr, degre
 from Analytic import Analytic_delay_vs_theta
@@ -23,50 +23,35 @@ def timing(time,weight,nbBins=100):
    binSize=dt[1:nbBins+1]-dt[0:nbBins]
    dNdt=dN/binSize *timecenter
-   return timecenter, dNdt
+   return timecenter, dNdt#/max(dNdt)
-def drawTiming(files,plot_case=0):
+def drawTiming(files,plot_gen=False,plot_energy_range=False):
    '''
       Plot flux versus time delay
       Input: list of directories
-      plot_case: (0) nothing, (1) all source spectrum, (2) all jets, (3) by energy band, 
-         (4) by generation 
       Output: graph of flux versus time delay 
    '''
    fig = figure()
    ax = fig.add_subplot(111)
    for fileId in files:
-      if plot_case == 1:
-         i=1
-         for powerlaw_index in [1,1.5,2,2.5]: 
-            delta_t,dNdt = ReadTiming(fileId,[0,i])
-            ax.plot(delta_t,dNdt,linestyle="steps-mid",label="p=%.1f"%powerlaw_index)
-            i+=1
-      elif plot_case == 2:
-         i=5
-         for jet_opening in ["isotrop","60 degre","30 degre"]: 
-            delta_t,dNdt = ReadTiming(fileId,[0,i])
-            ax.plot(delta_t,dNdt,linestyle="steps-mid",label=jet_opening)
-            i+=1
-      elif plot_case == 3:
-         i=8
-         Emin = [1e-3,1e0,1e3] #GeV
-         Emax = [1e0,1e3,1e5] #GeV
+      delta_t,dNdt = ReadTiming(fileId,[0,1])
+      p=ax.plot(delta_t,dNdt,linestyle="steps-mid")
+
+      if plot_energy_range:
+         labels = ["MeV band","GeV band","TeV band"]
+         i = 2
          for n in arange(0,3,1):
-            label = "%1.0e Gev"%Emin[n]+"< $E_\\gamma$ < %1.0e Gev"%Emax[n]
             delta_t,dNdt = ReadTiming(fileId,[0,i])
-            ax.plot(delta_t,dNdt,linestyle="steps-mid",label=label)
+            ax.plot(delta_t,dNdt,linestyle="steps-mid",label=labels[n])
             i+=1
-      elif plot_case == 4:
-         i=11
-         generation = ReadExtraFile(fileId,[3])
-         for gen in list(set(generation)):
+
+      if plot_gen:
+         generation = ReadGeneration(fileId,[0])
+         i=5
+         for gen in generation:
             delta_t,dNdt = ReadTiming(fileId,[0,i])
             ax.plot(delta_t,dNdt,linestyle="steps-mid",label="gen=%.0f"%gen)
             i+=1
-      else:
-         delta_t,dNdt = ReadTiming(fileId,[0,3])
-         ax.plot(delta_t,dNdt,linestyle="steps-mid",label=fileId)
    ax.set_xscale('log')
    ax.set_yscale('log')
@@ -81,12 +66,19 @@ def drawConvolveTiming(files):
    fig = figure()
    ax = fig.add_subplot(111)
    for fileId in files:
-      delta_t,dNdt1 = ReadTiming(fileId,[0,8])  # MeV band
-      delta_t,dNdt2 = ReadTiming(fileId,[0,9])  # GeV band
-      delta_t,dNdt3 = ReadTiming(fileId,[0,10]) # TeV band
-      ax.plot(delta_t,convolve(dNdt1,dNdt2,'same'),linestyle="steps-mid",label="MeV - GeV band")
-      ax.plot(delta_t,convolve(dNdt2,dNdt3,'same'),linestyle="steps-mid",label="GeV - TeV band")
-      ax.plot(delta_t,convolve(dNdt1,dNdt3,'same'),linestyle="steps-mid",label="MeV - TeV band")
+      delta_t,dNdt1 = ReadTiming(fileId,[0,2]) # MeV band
+      delta_t,dNdt2 = ReadTiming(fileId,[0,3]) # GeV band
+      delta_t,dNdt3 = ReadTiming(fileId,[0,4]) # TeV band
+      conv1 = convolve(dNdt1,dNdt2,'same')
+      conv2 = convolve(dNdt2,dNdt3,'same')
+      conv3 = convolve(dNdt1,dNdt3,'same')
+      ax.plot(delta_t,conv1/max(conv1),linestyle="steps-mid",label="MeV - GeV band")
+      ax.plot(delta_t,conv2/max(conv2),linestyle="steps-mid",label="GeV - TeV band")
+      ax.plot(delta_t,conv3/max(conv3),linestyle="steps-mid",label="MeV - TeV band")
+      print fileId, " ====================================="
+      print "   MeV - GeV band: ",delta_t[r_[True, conv1[1:] > conv1[:-1]] & r_[conv1[:-1] > conv1[1:], True]] 
+      print "   GeV - TeV band: ",delta_t[r_[True, conv2[1:] > conv2[:-1]] & r_[conv2[:-1] > conv2[1:], True]]
+      print "   MeV - TeV band: ",delta_t[r_[True, conv3[1:] > conv3[:-1]] & r_[conv3[:-1] > conv3[1:], True]]
    ax.set_xscale('log')
    ax.set_yscale('log')
@@ -117,21 +109,26 @@ def drawDelay_vs_angle(fileId):
    ax = fig.add_subplot(111)
    nbBins = 100
-   theta  = ReadExtraFile(fileId,[4])*degre
+   generation, theta = ReadExtraFile(fileId,cols=[3,4])
+   theta *= degre
    energy = ReadEnergy(fileId)
    delay  = ReadTime(fileId)
    Emin = [1e-3,1e0,1e3] #GeV
    Emax = [1e0,1e3,1e5] #GeV
-   for n in arange(0,3,1):
-      cond= (energy>Emin[n]) & (energy<Emax[n])
-      if Emin[n]==1e-3:
-         label="MeV band"
-      if Emin[n]==1e0:
-         label="GeV band"
-      if Emin[n]==1e3:
-         label="TeV band"
-      ax.plot(theta[cond],delay[cond],".",label=label)
+
+   for gen in list(set(generation)):
+      cond = generation==gen
+      ax.plot(theta[cond],delay[cond],",",label="gen = %.0f"%gen)
+   #for n in arange(0,3,1):
+   #   cond= (energy>Emin[n]) & (energy<Emax[n])
+   #   if Emin[n]==1e-3:
+   #      label="MeV band"
+   #   if Emin[n]==1e0:
+   #      label="GeV band"
+   #   if Emin[n]==1e3:
+   #      label="TeV band"
+   #   ax.plot(theta[cond],delay[cond],".",label=label)
    dt,angle=delay_vs_theta(theta,delay,nbBins)
    ax.plot(angle,dt,color='m',linestyle="steps-mid",linewidth=2)
@@ -0,0 +1 @@
+from Constants import *
@@ -0,0 +1,2 @@
+*
+!.gitignore
@@ -0,0 +1,57 @@
+#!/usr/bin/python
+
+from sys import argv
+import matplotlib.pyplot as plt
+from matplotlib import gridspec
+from numpy import *
+from scipy.integrate import quad
+from Modules.Analytic import Ethreshold_gg
+
+c=2.99792458e10 # cm.s-1
+Mpc=(3.0856776e+16)*1e8 # Mpc to cm
+H0=67.8*1e5/(Mpc) # s-1
+omegaM = 0.3
+omegaK = 0
+omegaL = 0.7
+zlim=-0.
+def properIntegrand(z):
+   return -c/(H0*(1+z)*sqrt(omegaM*(1+z)**3+omegaK*(1+z)**2+omegaL))
+def comobileIntegrand(z):
+   return -c/(H0*sqrt(omegaM*(1+z)**3+omegaK*(1+z)**2+omegaL))
+
+color=['b','r','g','m','c','y']
+
+# Proper distance figure
+#========================
+fig1 = plt.figure()
+ax11 = fig1.add_subplot(111)
+
+labels=["0","0.5","1","2","3"]
+theory=[1,2,3,4,5]
+ind=0
+for lab in labels:
+   # theoritical curve ===================================================
+   # f: without cosmo
+   # g: with cosmo
+   e,f,g=loadtxt('lambda_e.dat',unpack=True,usecols=[0,int(theory[ind]),int(theory[ind])+5])
+   e = e*511.e3/1.e9 #GeV
+   cond= (e>100) & (e<1e5)
+   e=e[cond]
+   g=g[cond]
+   f=f[cond]
+   p = ax11.plot(e,g,"-"+color[ind],label="z = "+lab) 
+   ax11.plot(e,f,color=p[0].get_color(),linestyle='--')
+
+   ind=ind+1
+
+ax11.axvline(x=Ethreshold_gg(), ymin=0., ymax = 1., color='r', linewidth=2)
+
+ax11.legend(loc="best")
+ax11.grid(b=True,which='major')
+ax11.set_ylim([1e-4,1e4])
+ax11.set_xscale('log')
+ax11.set_yscale('log')
+ax11.set_xlabel("energy [GeV]")
+ax11.set_ylabel("$\lambda_{\gamma\gamma}$ [Mpc]")
+
+plt.show()