arrival_angle.py
2.68 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
from numpy import log10, histogram, arange
from matplotlib.pyplot import figure, show
from matplotlib import gridspec
from read import ReadRadial, ReadAngleVsEnergy, ReadGeneration
from constants import degre
def arrivalAngle(theta,weight,nbBins=50,theta_range=[1e-6,180]):
cond = theta>0
theta = log10(theta[cond])
weight = weight[cond]
dn,angle2 = histogram(theta,nbBins,range=log10(theta_range),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 drawArrivalAngle(files,plot_gen=False,plot_energy_range=False,plot_time_range=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:
theta,dndtheta = ReadRadial(fileId,[0,1])
ax.plot(theta,theta*dndtheta,drawstyle='steps-mid',label=fileId)
if plot_energy_range:
theta,dNdtheta1,dNdtheta2,dNdtheta3 = ReadRadial(fileId,[0,2,3,4])
ax.plot(theta,theta*dNdtheta1,drawstyle='steps-mid',linestyle='--',label="MeV band")
ax.plot(theta,theta*dNdtheta2,drawstyle='steps-mid',linestyle='--',label="GeV band")
ax.plot(theta,theta*dNdtheta3,drawstyle='steps-mid',linestyle='--',label="TeV band")
if plot_time_range:
theta,dNdtheta1,dNdtheta2,dNdtheta3 = ReadRadial(fileId,[0,5,6,7])
ax.plot(theta,theta*dNdtheta1,drawstyle='steps-mid',linestyle='--',label="$t_{obs}$=30days")
ax.plot(theta,theta*dNdtheta2,drawstyle='steps-mid',linestyle='--',label="$t_{obs}$=1yr")
ax.plot(theta,theta*dNdtheta3,drawstyle='steps-mid',linestyle='--',label="$t_{obs}$=5yrs")
if plot_gen:
theta,dNdtheta1,dNdtheta2,dNdtheta3 = ReadRadial(fileId,[0,8,9,10])
ax.plot(theta,theta*dNdtheta1,drawstyle='steps-mid',linestyle='--',label="gen=1")
ax.plot(theta,theta*dNdtheta2,drawstyle='steps-mid',linestyle='--',label="gen=2")
ax.plot(theta,theta*dNdtheta3,drawstyle='steps-mid',linestyle='--',label="gen=3")
#i=8
#for gen in ReadGeneration(fileId,[0]):
# theta,dndtheta2 = ReadRadial(fileId,[0,i])
# ax.plot(theta,dndtheta2,drawstyle='steps-mid',linestyle='--',label="gen=%.0f"%gen)
# i+=1
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("$\\theta*dN/d\\theta$ [arbitrary]")
show()