Add direct printing of the maps in Analysis.py

Thomas Fitoussi
1 parent 6633967c
Showing 7 changed files with 91 additions and 86 deletions Show diff stats
Analysis.py
EBL-spectrums.py
Modules/Map.py
Modules/Map.pyc
Modules/Read.pyc
Modules/Spectrum.py
Modules/Spectrum.pyc
@@ -3,12 +3,15 @@
 from sys import argv
 from numpy import append, savetxt, shape, array, newaxis, zeros, arange
 from Modules.Read import ReadEnergy, ReadTime, ReadExtraFile, ReadProfile
+from Modules.Read import ReadMomentumAngle, ReadPositionAngle, resultDirectory
 from Modules.Spectrum import spectrum
+from Modules.Map import printMap, isotrop
 from Modules.Angle import radial, angle_vs_energy
 from Modules.Timing import timing
+from Modules.Constants import degre
-def InProgress(fileId,i,Nmax):
-   print "Work on ", fileId, ": ", (i*100)/Nmax, "% done"
+def InProgress(fileId,step,i,Nmax):
+   print "Work on ", fileId, "(step "+ step +": ", (i*100)/Nmax, "% done"
 if shape(argv)[0] < 2: 
    print "not enough arguments (at least 1)"
@@ -17,8 +20,10 @@ if shape(argv)[0] &lt; 2:
 #================================================================#
 #  Parameters
 #================================================================#
-folder = "Results/"
 PowerSpectrum=[1,1.5,2,2.5] 
+powerlaw_index = 2
+Elim = 1e-1 # GeV
+thetalim = 20 # degre
 for fileId in argv[1:]:
    # read files
@@ -27,16 +32,24 @@ for fileId in argv[1:]:
    weightini, generation, theta_arrival, Esource = ReadExtraFile(fileId,[2,3,4,5])
    nbPhotonsEmitted=ReadProfile(fileId,[3]) 
-   lim =100000
-   time=time[:lim]
-   energy=energy[:lim]
-   weightini=weightini[:lim]
-   generation = generation[:lim]
-   theta_arrival = theta_arrival[:lim]
-   Esource=Esource[:lim]
-
-   cond = (generation==0)
-   test=energy[cond]/Esource[cond]
+   #=============================================================================#
+   #  IMAGING 
+   #=============================================================================#
+   thetaDir,phiDir = ReadMomentumAngle(fileId)*degre
+   thetaPos,phiPos = ReadPositionAngle(fileId)*degre
+   theta = thetaDir - thetaPos
+   phi = phiDir -phiPos
+
+   # apply source spectrum
+   weight_source = (Esource/min(Esource))**(1-powerlaw_index)
+   weight = weightini* weight_source
+
+   # apply selection
+   cond=(energy>Elim) & (abs(theta)<thetalim) & (abs(phi)<thetalim)
+   weight=weight[cond]
+   theta=theta[cond]
+   phi=phi[cond]
+   printMap(theta,phi,weight,fileId,source=isotrop(),borne=[thetalim,thetalim])
    Gen_contrib = []
    Gen_cont = zeros((int(max(generation))+1))
@@ -48,14 +61,9 @@ for fileId in argv[1:]:
    for powerlaw_index in PowerSpectrum:
       # apply source spectrum
-      weight_source = Esource**(1-powerlaw_index)
+      weight_source = (Esource/min(Esource))**(1-powerlaw_index)
       weight = weightini* weight_source
-      print "%e" %(max(weight_source)/min(weight_source))
-
-      #cond = (energy>95) & (energy<105) & (generation==0)
-      #print weightini[cond], weight[cond]
-
       #=============================================================================#
       # GENERATION CONTRIBUTION
       #=============================================================================#
@@ -73,8 +81,7 @@ for fileId in argv[1:]:
       nbBins = 100
       # draw source spectrum
       Es=array(list(set(Esource)))
-      print shape(Es)
-      Ws= (Es)**(1-powerlaw_index)
+      Ws= (Es/min(Es))**(1-powerlaw_index)
       Es,Fs = spectrum(Es,Ws,nbBins=nbBins)
       Es=Es[:,newaxis]
       Fs=Fs[:,newaxis]
@@ -132,11 +139,11 @@ for fileId in argv[1:]:
       else:
          Timing = append(Timing,dNdt,axis=1)
-      InProgress(fileId,PowerSpectrum.index(powerlaw_index)+1,shape(PowerSpectrum)[0])
+      InProgress(fileId,"1",PowerSpectrum.index(powerlaw_index)+1,shape(PowerSpectrum)[0])
-   savetxt(folder+fileId+"/Generation.txt",Gen_contrib)
-   savetxt(folder+fileId+"/Spectrum.txt",Spectrum)
-   savetxt(folder+fileId+"/Source_spectrum.txt",Source)
-   savetxt(folder+fileId+"/Angle_vs_Energy.txt",Angle_Energy)
-   savetxt(folder+fileId+"/Radial_distribution.txt",Radial)
-   savetxt(folder+fileId+"/Timing.txt",Timing)
+   savetxt(resultDirectory+fileId+"/Generation.txt",Gen_contrib)
+   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)
@@ -5,8 +5,8 @@ from scipy.optimize import curve_fit
 import matplotlib.pyplot as plt
 from Modules.Constants import *
-fig = plt.figure()
-ax = fig.add_subplot(111)
+fig1 = plt.figure()
+ax = fig1.add_subplot(111)
 z=2
@@ -16,11 +16,11 @@ lamb,lambdaI = np.loadtxt(&quot;EBL_files/lambdaI_Dominguez.dat&quot;,unpack=True,usecols=
 hv = h*c/(lamb*1e-4)  # erg
 density = 1e-6*(4*pi/c)*lambdaI/(hv**2) /(erg_to_GeV*1e9) *(1+z)**3
 hv = hv *(erg_to_GeV*1e9)
-ax.plot(hv,density*hv**2,"--b",label="Dominguez and Al")
+ax.plot(hv,density*hv**2,"--b",label="Dominguez et Al")
 # ==== Kneiske and Doll - "best fit" ====
 hv,density = np.loadtxt("EBL_files/n_bestfit10.dat",unpack=True,usecols=[0,1])
-ax.plot(hv,density*hv**2,"--r",label="Kneiske and Doll - 'best fit'")
+ax.plot(hv,density*hv**2,"--r",label="Kneiske et Doll - 'best fit'")
 # ==== Kneiske and Doll - "lower limit" ====
 hv,density = np.loadtxt("EBL_files/n_lowerlimit10.dat",unpack=True,usecols=[0,179])
@@ -38,28 +38,50 @@ ax.plot(hv,density*hv**2,&quot;--c&quot;,label=&quot;Fraceschini&quot;)
 hv,density = np.loadtxt("EBL_files/n_Finke.dat",unpack=True,usecols=[0,201])
 #hv,density = np.loadtxt("EBL_files/n_Finke.dat",unpack=True,usecols=[0,1])
 density=density*(1+z)**3
-ax.plot(hv,density*hv**2,"--m",label="Finke and Al")
+ax.plot(hv,density*hv**2,"--m",label="Finke et Al")
 # ==== Gilmore ====
 hv,density = np.loadtxt("EBL_files/n_Gil.dat",unpack=True,usecols=[0,14])
 #hv,density = np.loadtxt("EBL_files/n_Gil.dat",unpack=True,usecols=[0,1])
-ax.plot(hv,density*hv**2,"--y",label="Gilmore and Al")
+ax.plot(hv,density*hv**2,"--y",label="Gilmore et Al")
 #==== CMB ====
-def nCMB(E):
+def nCMB(E,z):
    kTcmb = k*Tcmb*erg_to_GeV*1e9*(1+z)
    theta = E/kTcmb
    nCMB=(hb*c*erg_to_GeV*1e9)**(-3) *(E/np.pi)**2 /(np.exp(theta)-1)    
    return nCMB
 hv = np.logspace(-4,-1,1000)
-ax.plot(hv,nCMB(hv)*hv**2,"-k",label="CMB")
+ax.plot(hv,nCMB(hv,z)*hv**2,"-k",label="CMB")
 ax.set_xscale('log')
 ax.set_yscale('log')
 ax.grid(b=True,which='major')
-ax.legend(loc="best",title="z = %.0f"%z)
+ax.legend(loc="best",frameon=False,framealpha=0.5)
 ax.set_xlabel("energy [eV]")
-ax.set_ylabel("$n_{EBL}$ [photon.eV.cm$^{-3}$]")
+ax.set_ylabel("$n$ [photon.eV.cm$^{-3}$]")
+
+#==== Dominguez vs redshift ==================================================================
+fig2 = plt.figure()
+ax2 = fig2.add_subplot(111)
+
+z = np.loadtxt("EBL_files/z_Dominguez.dat",unpack=True,usecols=[0])
+for z_index in [1,9,12,15,17]:
+   lamb,lambdaI = np.loadtxt("EBL_files/lambdaI_Dominguez.dat",unpack=True,usecols=[0,z_index])
+   hv = h*c/(lamb*1e-4)  # erg
+   density = 1e-6*(4*pi/c)*lambdaI/(hv**2) /(erg_to_GeV*1e9) *(1+z[z_index])**3
+   hv = hv *(erg_to_GeV*1e9)
+   p=ax2.plot(hv,density*hv**2,"--")
+   # CMB
+   hv = np.logspace(-4,-1,1000)
+   ax2.plot(hv,nCMB(hv,z[z_index])*hv**2,color=p[0].get_color(),linestyle='-',label="z="+str(z[z_index-1]))
+
+ax2.set_xscale('log')
+ax2.set_yscale('log')
+ax2.grid(b=True,which='major')
+ax2.legend(loc="best",frameon=False,framealpha=0.5)
+ax2.set_xlabel("energy [eV]")
+ax2.set_ylabel("$n$ [photon.eV.cm$^{-3}$]")
 plt.show()
 from numpy import shape, pi, histogram2d, arange, newaxis, exp, log, array, sin, cos
 from scipy.signal import convolve2d
-from matplotlib.pyplot import figure, show, matshow
+from matplotlib.pyplot import figure, show, matshow, savefig
 from matplotlib.colors import LogNorm
-from Read import ReadEnergy, ReadExtraFile, ReadMomentumAngle, ReadPositionAngle, ReadProfile
+from Read import  resultDirectory
 from Constants import degre
 thetalim = 180
@@ -14,15 +14,15 @@ thetalim=float(thetalim)
 limits = [[-thetalim,thetalim],[-philim,philim]]
 def isotrop():
-   return array([[1]])
+   return array([[1]]), "isotropic"
 def Gaussian2D(center=[0,0],fwhm=0.01):
    x = arange(-thetalim,thetalim+step,step)
    y = arange(-philim,philim+step,step)
    y = y[:,newaxis]
-   return exp(-4*log(2) * ((x-center[0])**2 + (y+center[1])**2) / fwhm**2)
+   return exp(-4*log(2) * ((x-center[0])**2 + (y+center[1])**2) / fwhm**2), "jet_"+str(fwhm)
-def drawMap(files,Elim=0.1,source=isotrop(),borne=[180,180]):
+def printMap(theta,phi,weight,fileId,source=isotrop(),borne=[180,180]):
    '''
       Plot 2D histogram of arrival direction of the photons with energy upper than Elim
       Input: list of directories
@@ -31,47 +31,23 @@ def drawMap(files,Elim=0.1,source=isotrop(),borne=[180,180]):
       Input (optional): limits on theta, phi 
       Output: 2D histogram of theta cos(phi), theta sin(phi)
    '''
-   print "Cut at", Elim, "GeV" 
    fig = figure()
-   nbplots = len(files)
-   ind=1
-   for fileId in files:
-      thetaDir,phiDir = ReadMomentumAngle(fileId)*degre
-      thetaPos,phiPos = ReadPositionAngle(fileId)*degre
-      theta = thetaDir - thetaPos
-      phi = phiDir -phiPos
-
-      energy = ReadEnergy(fileId)
-      weightini, Esource = ReadExtraFile(fileId,[2,5])
-      
-      zsource, Dsource, nbPhotonsEmitted, EGMF = ReadProfile(fileId,[1,2,3,4])
-      # apply source spectrum
-      powerlaw_index = 2
-      weight = weightini*(Esource/min(Esource))**(1-powerlaw_index)
-
-      cond=(energy>Elim) & (abs(theta)<thetalim) & (abs(phi)<philim)
-      energy = energy[cond]
-      weight=weight[cond]/nbPhotonsEmitted
-      theta=theta[cond]
-      phi=phi[cond]
-
-      H,xedges,yedges = histogram2d(theta*cos(phi),theta*sin(phi),nbBins,weights=weight,range=limits)
-      conv=convolve2d(H,source,boundary="wrap",mode="same")
-
-      ax = fig.add_subplot(100+nbplots*10+ind)
-      im = ax.matshow(conv,norm=LogNorm(),aspect='auto',
-                        extent=[xedges[0],xedges[-1],yedges[0],yedges[-1]])
-      #H,xedges,yedges,im = ax.hist2d(theta,phi,nbBins,weights=weight,norm=LogNorm())
-      cbar=fig.colorbar(im, ax=ax)
-      cbar.ax.set_ylabel("counts normalize to 1 photon emitted")
-      ax.set_xlim((-borne[0],borne[0]))
-      ax.set_ylim((-borne[1],borne[1]))
-      ax.set_xlabel("$\\theta$ $\\cos(\\phi)$ [deg]")
-      ax.set_ylabel("$\\theta$ $\\sin(\\phi)$ [deg]")
-      ax.grid(b=True,which='major')
-      ax.set_title("source at %.0f"%Dsource+"Mpc (z=%.4f"%zsource+"): EGMF="+str(EGMF)+" Gauss"   )
-    
-      ind=ind+1
-
-   show()
+   H,xedges,yedges = histogram2d(theta*cos(phi),theta*sin(phi),nbBins,weights=weight,range=limits)
+   conv=convolve2d(H,source[0],boundary="wrap",mode="same")
+
+   ax = fig.add_subplot(111)
+   im = ax.matshow(conv,norm=LogNorm(),aspect='auto',
+                     extent=[xedges[0],xedges[-1],yedges[0],yedges[-1]])
+   #H,xedges,yedges,im = ax.hist2d(theta,phi,nbBins,weights=weight,norm=LogNorm())
+   cbar=fig.colorbar(im, ax=ax)
+   cbar.ax.set_ylabel("counts")
+   ax.set_xlim((-borne[0],borne[0]))
+   ax.set_ylim((-borne[1],borne[1]))
+   ax.set_xlabel("$\\theta$ $\\cos(\\phi)$ [deg]")
+   ax.set_ylabel("$\\theta$ $\\sin(\\phi)$ [deg]")
+   ax.grid(b=True,which='major')
+   ax.set_title(fileId+" - "+source[1])
+
+   savefig(resultDirectory+fileId+"/map_"+source[1]+".png")
+   #show()
@@ -51,7 +51,7 @@ def drawSpectrum(files,PlotAnalytic=False):
       for powerlaw_index in [1,1.5,2,2.5]: 
          # read files
          Es,Fs = ReadSourceSpectrum(fileId,[0,i])
-         p = ax1.plot(Es,Fs,drawstyle='.')
+         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)