Angle.py
4.28 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
131
132
133
134
135
#!/usr/bin/python
from sys import argv
import matplotlib.pyplot as plt
from matplotlib import gridspec
from Constantes import *
if argv[1] == "EGMF":
fileName="Results_EGMF"
elif argv[1] == "EBL":
fileName="Results_EBL"
B=10**(-15)
Model=["Kneiske and Doll - 'best fit'","Kneiske and Doll - 'lower limit'",
"Fraceschini","Finke and Al","Gilmore and Al","Dominguez and Al"]
else:
print "Used: ./Angle_distribution.py arg1 ind1 ind2 ..."
print " arg1 = EGMF or EBL "
print " ind1 = file number \n"
print "Examples: "
print " to study EGMF14: ./Angle_distribution.py EGMF 14 "
print " to study EBL1 and EBL2: ./Angle_distribution.py EBL 1 2 \n"
exit()
nbBins = 100
Eliminf = 1e0 #GeV
Elimsup = 1e5 #GeV
Esource=100 #TeV
Dsource=135 #Mpc
delta_to_theta=(lambda_gg/Esource*1e3)/Dsource
def fit_theta(E,B):
RL0=RL*(Esource/2)*(1e-17/B)
Dic0=Dic/(Esource/2)
delta=Dic0/(2*RL0)*((Esource/2)**2 *Eic/E-1)
return abs(arcsin(delta_to_theta*sin(delta))*degre)
fig2 = plt.figure()
ax2 = fig2.add_subplot(111)
fig1 = plt.figure()
gs = gridspec.GridSpec(2, 1, height_ratios=[4,1])
fig1.subplots_adjust(hspace=0.001)
ax11 = fig1.add_subplot(gs[0])
ax12 = fig1.add_subplot(gs[1],sharex=ax11)
color=['b','r','g','c','m','y']
ind = 0
for fileId in argv[2:]:
energy,dirx,diry,dirz = loadtxt(fileName+fileId+"/Results_momentum",unpack=True,usecols=[0,1,2,3])
posx,posy,posz = loadtxt(fileName+fileId+"/Results_position",unpack=True,usecols=[1,2,3])
charge,weight,gen=loadtxt(fileName+fileId+"/Results_extra",unpack=True,usecols=[0,2,3])
hyp=sqrt(posx**2+posy**2+posz**2)
posx=posx/hyp
posy=posy/hyp
posz=posz/hyp
hyp=sqrt(dirx**2+diry**2+dirz**2)
dirx=dirx/hyp
diry=diry/hyp
dirz=dirz/hyp
costheta=dirx*posx+diry*posy+dirz*posz
cond=(costheta<=1)&(costheta>=-1)&(energy>Eliminf)&(energy<Elimsup)
theta = arccos(costheta[cond])*degre
weight= weight[cond]
energy= energy[cond]
print fileId, ":", shape(energy)
# Angle versus energy
#=====================
nbBins = 40
Edata =log10(energy)
angle,ener =histogram(Edata,nbBins,weights=theta)
nb,ener =histogram(Edata,nbBins)
ener = 10**ener
enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2
angle = angle/nb
if argv[1] == "EGMF":
ax11.plot(enercenter,angle,'.'+color[ind],
label="$10^{-%.0f"%float(fileId)+"}$Gauss - MC")
#ax11.plot(energy,theta,'.'+color[ind],label="$10^{-%.0f"%float(fileId)+"}$Gauss - MC")
B=10**(-float(fileId))
yfit = fit_theta(enercenter,B)
ax11.plot(enercenter,yfit,'--'+color[ind],linewidth=2,label="$10^{-%.0f"%float(fileId)+"}$Gauss - analytic")
elif argv[1] == "EBL":
ax11.plot(enercenter,angle,'.'+color[ind],label=Model[int(fileId)-1])
if ind==0:
yfit = fit_theta(enercenter,B)
ax11.plot(enercenter,yfit,'--'+color[ind],linewidth=2,label="$10^{-15}$Gauss - analytic")
error=angle/yfit-1
ax12.plot(enercenter,error,'+'+color[ind])
# figure: radial distribution
#=============================
nbBins = 20
#theta = theta[theta!=0]
#weight=weight[theta!=0]
#theta = log10(theta**2)
dn,angle1=histogram(theta**2,nbBins, range=[0,0.25],weights=weight)
#angle1=10**angle1
thetacenter=(angle1[1:nbBins+1]+angle1[0:nbBins])/2
dtheta=angle1[1:nbBins+1]-angle1[0:nbBins]
dndtheta2=dn/dtheta
ax2.plot(thetacenter,dndtheta2,'-'+color[ind],
drawstyle='steps-mid',label="$10^{-%.0f"%float(fileId)+"}$Gauss - MC")
#ax2.hist(theta, nbBins, range=[0,0.5], weights=weight,log=True,
# facecolor=color[ind],alpha=.5,label="$10^{-%.0f"%float(fileId)+"}$Gauss")
ind=ind+1
ax11.legend(loc="best")
ax11.set_xscale('log')
ax11.set_yscale('log')
ax11.set_ylim([0,180])
ax11.grid(b=True,which='major')
ax11.set_ylabel("$\\theta$ [degre]")
ax12.set_xscale('log')
ax12.grid(b=True,which='major')
ax12.set_xlabel("energy [GeV]")
ax12.set_ylabel("relative error")
step=0.2*(max(error)-min(error))
plt.ylim(min(error)-step, max(error)+step)
xticklabels = ax11.get_xticklabels()
plt.setp(xticklabels, visible=False)
ax2.legend(loc="best",title="%3d - %3d GeV"%(Eliminf,Elimsup))
#ax2.set_xscale('log')
ax2.set_yscale('log')
ax2.grid(b=True,which='major')
ax2.set_xlabel("$\\theta^2$ [deg$^2$]")
ax2.set_ylabel("$dN/d\\theta^2$ [deg$^{-2}$]")
plt.show()