Map: convolution of the results from mono-directional source

-> source isotrop -> gaussian

Map: convolution of the results from mono-directional source
-> source isotrop -> gaussian
Thomas Fitoussi
1 parent 374f3d52
Showing 9 changed files with 81 additions and 62 deletions Show diff stats
Analysis.py
Modules/Angle.py
Modules/Angle.pyc
Modules/Map.py
Modules/Map.pyc
Modules/Spectrum.py
Modules/Spectrum.pyc
Modules/Timing.py
Modules/Timing.pyc
@@ -6,7 +6,7 @@ from Modules.Spectrum import drawSpectrum, drawSpectrum_normaLepton
 from Modules.Angle import drawAngle, drawRadial
 from Modules.Timing import drawTiming
 from Modules.Generation import drawGeneration
-from Modules.Map import drawMap, drawMapHisto
+from Modules.Map import drawMap
  
 def error(error_type):
    print error_type
@@ -43,9 +43,6 @@ elif argv[1] == &quot;timing&quot;:
 elif argv[1] == "generation":
    drawGeneration(argv[2:])
  
-elif argv[1] == "maphisto":
-   drawMapHisto(argv[2:])
-
 elif argv[1] == "map":
    drawMap(argv[2:])
  
@@ -83,6 +83,8 @@ def radial(fileName):
    thetacenter=(angle[1:nbBins+1]+angle[0:nbBins])/2
    dtheta=angle[1:nbBins+1]-angle[0:nbBins]
    dndtheta2=dn/dtheta
+   maxflux = ReadNphot_source(fileName)
+   dndtheta2 /= maxflux
  
    return thetacenter, dndtheta2
  
@@ -92,8 +94,7 @@ def drawRadial(files):
  
    for fileId in files:
       theta,dndtheta2=radial(fileId)
-      Nmax=ReadNphot_source(fileId) # Nb photons emitted by the source
-      ax.plot(theta,dndtheta2/Nmax,drawstyle='steps-mid',label=fileId)
+      ax.plot(theta,dndtheta2,drawstyle='steps-mid',label=fileId)
  
    ax.legend(loc="best",title="%3d - %3d GeV"%(Eliminf,Elimsup))
    ax.set_yscale('log')
-from numpy import shape, pi
-from matplotlib.pyplot import figure, show
+from numpy import shape, pi, histogram2d, arange, newaxis, exp, log, array
+from scipy.signal import convolve2d
+from matplotlib.pyplot import figure, show, matshow
 from matplotlib.colors import LogNorm
 from Read import ReadEnergy, ReadWeight, ReadMomentumAngle, ReadPositionAngle, ReadNphot_source
 from Read import ReadD_source, Readz_source, ReadEGMF
 from Constants import degre
  
-Elim=1e-1 # GeV
-borne=5 # degree
-nbBins = 100
+thetalim = 180
+philim = 90
+step = 1 #deg
+nbBins = array([int(2*philim/step)+1,int(2*thetalim/step)+1])
+philim=float(philim)
+thetalim=float(thetalim)
+limits = [[-thetalim,thetalim],[-philim,philim]]
  
-def thetaphi(fileId):
+def thetaphi(fileId,Elim=0.1):
    energy = ReadEnergy(fileId)
    weight = ReadWeight(fileId)
    thetaDir,phiDir = ReadMomentumAngle(fileId)*degre
    thetaPos,phiPos = ReadPositionAngle(fileId)*degre
+   Nmax=ReadNphot_source(fileId) # Nb photons emitted by the source
  
    theta = thetaDir - thetaPos
    phi = phiDir -phiPos
  
-   cond=(energy>Elim) & (abs(theta)<borne) & (abs(phi)<borne)
+   cond=(energy>Elim) & (abs(theta)<thetalim) & (abs(phi)<philim)
    energy = energy[cond]
-   weight=weight[cond]
+   weight=weight[cond]/Nmax
    theta=theta[cond]
    phi=phi[cond]
  
    print "file", fileId, "->", shape(energy)[0], "events"
    return energy, theta, phi, weight
  
-def drawMap(files):
-   print "Cut at", Elim, "GeV" 
-   fig = figure()
-   nbplots = len(files)
+def isotrop():
+   return array([[1]])
  
-   ind=1
-   for fileId in files:
-      energy, theta, phi, weight = thetaphi(fileId)
+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)
  
-      ax = fig.add_subplot(100+nbplots*10+ind)
-
-      Emax = max(energy)
-      Emin = min(energy)
-      colors = energy/Emax
-      area = pi * (15 * colors)**2 +15
-      sax = ax.scatter(theta, phi, c=colors, s=area, alpha=0.3, marker='o')
-
-      if Emin < 1e-3:
-         liminf = "%.0f"%(Emin*1e6)+"keV"
-      elif Emin < 1:
-         liminf = "%.0f"%(Emin*1e3)+"MeV"
-      else: 
-         liminf = "%.0f"%(Emin)+"GeV"
-      limsup = "%.0f"%(Emax*1e-3)+"TeV"
-      cbar = fig.colorbar(sax, ticks=[0.01,1])
-      cbar.ax.set_yticklabels([liminf,limsup])
-
-      ax.set_xlim((-borne,borne))
-      ax.set_ylim((-borne,borne))
-      ax.set_xlabel("$\\theta$ [deg]")
-      ax.set_ylabel("$\\phi$ [deg]")
-      zsource=Readz_source(fileId)
-      Dsource=ReadD_source(fileId)
-      EGMF=ReadEGMF(fileId)
-      ax.set_title("source at %.0f"%Dsource+"Mpc (z=%.4f"%zsource+"): EGMF="+str(EGMF)+" Gauss"   )
- 
-      ind=ind+1
-
-   show()
-
-def drawMapHisto(files):
+def drawMap(files,Elim=0.1,source=isotrop(),borne=[180,90]):
    print "Cut at", Elim, "GeV" 
    fig = figure()
    nbplots = len(files)
  
    ind=1
    for fileId in files:
-      energy, theta, phi, weight = thetaphi(fileId)
-      Nmax=ReadNphot_source(fileId) # Nb photons emitted by the source
+      energy, theta, phi, weight = thetaphi(fileId,Elim)
+
+      H,xedges,yedges = histogram2d(phi,theta,nbBins,weights=weight,range=limits)
+      conv=convolve2d(H,source,boundary="wrap",mode="same")
  
       ax = fig.add_subplot(100+nbplots*10+ind)
-      H,xedges,yedges,im = ax.hist2d(theta,phi,nbBins,range=[[-borne,borne],[-borne,borne]],weights=weight/Nmax,norm=LogNorm())
-      extent = [yedges[0], yedges[-1], xedges[0], xedges[-1]]
+      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,borne))
-      ax.set_ylim((-borne,borne))
+      ax.set_xlim((-borne[0],borne[0]))
+      ax.set_ylim((-borne[1],borne[1]))
       ax.set_xlabel("$\\theta$ [deg]")
       ax.set_ylabel("$\\phi$ [deg]")
+      ax.grid(b=True,which='major')
       zsource=Readz_source(fileId)
       Dsource=ReadD_source(fileId)
       EGMF=ReadEGMF(fileId)
@@ -94,3 +73,43 @@ def drawMapHisto(files):
       ind=ind+1
  
    show()
+
+#def drawMap(files):
+#   print "Cut at", Elim, "GeV" 
+#   fig = figure()
+#   nbplots = len(files)
+#
+#   ind=1
+#   for fileId in files:
+#      energy, theta, phi, weight = thetaphi(fileId)
+#
+#      ax = fig.add_subplot(100+nbplots*10+ind)
+#
+#      Emax = max(energy)
+#      Emin = min(energy)
+#      colors = energy/Emax
+#      area = pi * (15 * colors)**2 +15
+#      sax = ax.scatter(theta, phi, c=colors, s=area, alpha=0.3, marker='o')
+#
+#      if Emin < 1e-3:
+#         liminf = "%.0f"%(Emin*1e6)+"keV"
+#      elif Emin < 1:
+#         liminf = "%.0f"%(Emin*1e3)+"MeV"
+#      else: 
+#         liminf = "%.0f"%(Emin)+"GeV"
+#      limsup = "%.0f"%(Emax*1e-3)+"TeV"
+#      cbar = fig.colorbar(sax, ticks=[0.01,1])
+#      cbar.ax.set_yticklabels([liminf,limsup])
+#
+#      ax.set_xlim((-borne,borne))
+#      ax.set_ylim((-borne,borne))
+#      ax.set_xlabel("$\\theta$ [deg]")
+#      ax.set_ylabel("$\\phi$ [deg]")
+#      zsource=Readz_source(fileId)
+#      Dsource=ReadD_source(fileId)
+#      EGMF=ReadEGMF(fileId)
+#      ax.set_title("source at %.0f"%Dsource+"Mpc (z=%.4f"%zsource+"): EGMF="+str(EGMF)+" Gauss"   )
+# 
+#      ind=ind+1
+#
+#   show()
@@ -47,7 +47,7 @@ def drawSpectrum(files,PlotElmag=False,PlotAnalytic=False):
          ax.plot(energy,flux,"--",color=p[0].get_color(),label="Elmag - "+fileId)
  
    if PlotAnalytic:
-      alpha=5e3
+      alpha=2.8e3
       ax.plot(energy,alpha*sqrt(energy),'--m',linewidth=2)
       ax.plot(energy,energy**(0)*max(flux),'--m',linewidth=2)
       ax.axvline(x=Ethreshold_gg(), ymin=0., ymax = 1., color='r', linewidth=2)
@@ -9,7 +9,7 @@ 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]), "~", prop
+   print "file", fileName, "->", shape(time)[0], "events", "negative time:", shape(time[time<0])[0], "~", prop
  
    time=time[time>0]
    weight=weight[time>0]
@@ -17,13 +17,15 @@ def timing(fileName):
    nbBins = 100
    zSource = Readz_source(fileName)
    tlim = quad(comobileTime,zSource,0)[0]/yr
-   time=log10(time/tlim)
+   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