Map.py
2.82 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
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, ReadExtraFile, ReadMomentumAngle, ReadPositionAngle, ReadProfile
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 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:
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()