Map.py
2.22 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
from numpy import shape, histogram2d, arange, newaxis, array, nditer, zeros, histogram
from numpy import sqrt, pi, exp, log, log10, sin, cos, where
from matplotlib.pyplot import figure, matshow, savefig, show
from matplotlib.colors import LogNorm
#from scipy.signal import convolve2d
from Read import resultDirectory, ReadRadial, ReadAngleVsEnergy
from Analytic import Analytic_theta
from Constants import degre
thetalim = 180
philim = 180
step = 0.25 #deg
bining = array([int(2*philim/step)+1,int(2*thetalim/step)+1])
philim=float(philim)
thetalim=float(thetalim)
limits = [[-thetalim,thetalim],[-philim,philim]]
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
Input (optional): lower energy limit (in GeV (default = 100MeV)
Input (optional): nature of the source (default isotrop)
Input (optional): limits on theta, phi
Output: 2D histogram of theta cos(phi), theta sin(phi)
'''
fig = figure()
#===============================================================================================#
# IMAGE
#===============================================================================================#
ax = fig.add_subplot(111)
#H,xedges,yedges = histogram2d(theta*cos(phi),theta*sin(phi),bining,weights=weight,range=limits)
#conv=convolve2d(H,source[0],boundary="wrap",mode="same")
H,xedges,yedges,im = ax.hist2d(theta*cos(phi),theta*sin(phi),bining,weights=weight,
norm=LogNorm(),range=limits)
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(saveId+" - isotrop")
savefig(resultDirectory+saveId+"/map_isotrop.png")
#def isotrop():
# 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), "jet_"+str(fwhm)