Angle.py
3.43 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
from numpy import log10, histogram
from matplotlib.pyplot import figure, show
from matplotlib import gridspec
from Read import ReadRadial, ReadAngleVsEnergy
from Analytic import Analytic_theta
from Constants import degre
def angle_vs_energy(theta,energy,weight,nbBins=50):
'''
Compute arrival angle versus energy
Input: directory name
Output: energy, arrival angle
'''
energy =log10(energy)
angle,ener = histogram(energy,nbBins,weights=weight*theta)
nb, ener = histogram(energy,nbBins,weights=weight)
ener = 10**ener
enercenter = (ener[1:nbBins+1]+ener[0:nbBins])/2
angle /= nb
return enercenter, angle
def radial(theta,weight,nbBins=50):
cond = theta!=0
theta = log10(theta[cond])
weight = weight[cond]
dn,angle2 = histogram(theta,nbBins,range=[log10(1e-3),log10(25)],weights=weight)
angle2=10**angle2
#theta = theta**2
#dn,angle2 = histogram(theta,nbBins,range=[0,25],weights=weight)
thetacenter = (angle2[1:nbBins+1]+angle2[0:nbBins])/2
dtheta = angle2[1:nbBins+1]-angle2[0:nbBins]
dndtheta = dn/dtheta
dndtheta /= max(dndtheta)
return thetacenter, dndtheta
def drawRadial(files,all_source_spectrum=False,all_jets=False):
'''
Plot flux versus arrival angle
Input: list of directories
Output: graph of flux versus angle
'''
fig = figure()
ax = fig.add_subplot(111)
for fileId in files:
if all_source_spectrum:
i=1
for powerlaw_index in [1,1.5,2,2.5]:
theta,dndtheta2 = ReadRadial(fileId,[0,i])
ax.plot(theta,dndtheta2,drawstyle='steps-mid',label="p=%.1f"%powerlaw_index)
i+=1
elif all_jets:
i=5
for jet_opening in ["isotrop","60 deg","30 deg"]:
theta,dndtheta2 = ReadRadial(fileId,[0,i])
ax.plot(theta,dndtheta2,drawstyle='steps-mid',label=jet_opening)
i+=1
else:
theta,dndtheta2 = ReadRadial(fileId,[0,5])
ax.plot(theta,dndtheta2,drawstyle='steps-mid',label=fileId)
ax.legend(loc="best")
ax.set_xscale('log')
ax.set_yscale('log')
ax.grid(b=True,which='major')
ax.set_xlabel("$\\theta$ [deg]")
ax.set_ylabel("$dN/d\\theta$ [arbitrary]")
show()
def drawAngle_vs_energy(files,all_source_spectrum=False,all_jets=False,PlotAnalytic=False):
'''
Plot arrival angle versus energy + analytic expression
Input: list of directories
Output: graph of arrival angle versus energy
'''
fig = figure()
ax1 = fig.add_subplot(111)
for fileId in files:
if all_source_spectrum:
i=1
for powerlaw_index in [1,1.5,2,2.5]:
ener,angle = ReadAngleVsEnergy(fileId,[0,i])
p=ax1.plot(ener,angle,drawstyle='steps-mid',label="p=%.1f"%powerlaw_index)
elif all_jets:
i=5
for jet_opening in ["isotrop","60 deg","30 deg"]:
ener,angle = ReadAngleVsEnergy(fileId,[0,i])
p=ax1.plot(ener,angle,drawstyle='steps-mid',label=jet_opening)
else:
ener,angle = ReadAngleVsEnergy(fileId,[0,1])
p=ax1.plot(ener,angle,drawstyle='steps-mid',label=fileId)
if PlotAnalytic:
yfit = Analytic_theta(ener,fileId)
ax1.plot(ener,yfit,'--',color=p[0].get_color(),linewidth=2)
ax1.legend(loc="best")
ax1.set_xscale('log')
ax1.set_yscale('log')
ax1.set_ylim([0,180])
ax1.grid(b=True,which='major')
ax1.set_ylabel("$\\theta$ [deg]")
ax1.set_xlabel("energy [GeV]")
show()