Angle.py
1.4 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
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
Eliminf = 1e0 #GeV
Elimsup = 1e5 #GeV
def angle_vs_energy(theta,energy,weight,nbPhotonsEmitted=1,nbBins=40):
'''
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 * nbPhotonsEmitted
return enercenter, angle
def drawAngle_vs_energy(files):
'''
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:
ener,angle = ReadAngleVsEnergy(fileId,[0,1])
p=ax1.plot(ener,angle,drawstyle='steps-mid',label=fileId)
yfit = Analytic_theta(ener,fileId)
ax1.plot(ener,yfit,'--',color="m",linewidth=2,label="analytic")
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$ [rad]")
ax1.set_xlabel("energy [GeV]")
show()