Map.py
4.34 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
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
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.colors import LogNorm
from Read import ReadEnergy, ReadWeight, ReadMomentumAngle, ReadPositionAngle, ReadNphot_source
from Read import ReadD_source, Readz_source, ReadEGMF
from Constants import degre
thetalim = 180
philim = 180
step = 0.25 #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,Elim=0.1):
'''
Compute arrival angles (theta,phi) selected with energy upper than Elim
normalize to 1 source photon
Input: directory name
Input (optional): lower energy limit in GeV (default = 100MeV)
Output: energy (GeV), theta, phi (degre), weight
'''
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)<thetalim) & (abs(phi)<philim)
energy = energy[cond]
weight=weight[cond]/Nmax
theta=theta[cond]
phi=phi[cond]
print "file", fileId, "->", shape(energy)[0], "events"
return energy, theta, phi, weight
def isotrop():
return array([[1]])
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)
def drawMap(files,Elim=0.1,source=isotrop(),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)
'''
print "Cut at", Elim, "GeV"
fig = figure()
nbplots = len(files)
ind=1
for fileId in files:
energy, theta, phi, weight = thetaphi(fileId,Elim)
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')
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 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()