Angle.py
3.05 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
from numpy import sqrt, arccos, shape, log10, histogram
from matplotlib.pyplot import figure, show, ylim, setp
from matplotlib import gridspec
from Read import ReadMomentum, ReadPosition, ReadEnergy, ReadWeight, ReadNphot_source
from Analytic import Analytic_theta
from Constants import degre
Eliminf = 1e0 #GeV
Elimsup = 1e5 #GeV
def compute_theta(fileName):
position = ReadPosition(fileName)
momentum = ReadMomentum(fileName)
energy = ReadEnergy(fileName)
weight = ReadWeight(fileName)
hyp=sqrt(position[0]**2+position[1]**2+position[2]**2)
position /= hyp
hyp=sqrt(momentum[0]**2+momentum[1]**2+momentum[2]**2)
momentum /= hyp
costheta=position[0]*momentum[0]+position[1]*momentum[1]+position[2]*momentum[2]
cond=(costheta<=1)&(costheta>=-1)&(energy>Eliminf)&(energy<Elimsup)
theta = arccos(costheta[cond])*degre
weight= weight[cond]
energy= energy[cond]
print "file", fileName, "->", shape(energy)[0], "events"
return energy, weight, theta
def angle_vs_energy(fileName):
energy,weight,theta = compute_theta(fileName)
nbBins = 40
energy =log10(energy)
angle,ener =histogram(energy,nbBins,weights=theta)
nb,ener =histogram(energy,nbBins)
ener = 10**ener
enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2
angle = angle/nb
return enercenter, angle
def drawAngle(files):
fig = figure()
gs = gridspec.GridSpec(2, 1, height_ratios=[4,1])
fig.subplots_adjust(hspace=0.001)
ax1 = fig.add_subplot(gs[0])
ax2 = fig.add_subplot(gs[1],sharex=ax1)
for fileId in files:
energy,angle = angle_vs_energy(fileId)
p=ax1.plot(energy,angle,'.',label=fileId)
yfit = Analytic_theta(energy,fileId)
ax1.plot(energy,yfit,'--',color=p[0].get_color(),linewidth=2,label="analytic")
error=angle/yfit-1
ax2.plot(energy,error,'+',color=p[0].get_color())
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$ [degre]")
ax2.set_xscale('log')
ax2.grid(b=True,which='major')
ax2.set_xlabel("energy [GeV]")
ax2.set_ylabel("error")
step=0.2*(max(error)-min(error))
ylim(min(error)-step, max(error)+step)
xticklabels = ax1.get_xticklabels()
setp(xticklabels, visible=False)
show()
def radial(fileName):
nbBins = 20
energy,weight,theta = compute_theta(fileName)
dn,angle=histogram(theta**2,nbBins,range=[0,0.25],weights=weight)
thetacenter=(angle[1:nbBins+1]+angle[0:nbBins])/2
dtheta=angle[1:nbBins+1]-angle[0:nbBins]
dndtheta2=dn/dtheta
maxflux = ReadNphot_source(fileName)
dndtheta2 /= maxflux
return thetacenter, dndtheta2
def drawRadial(files):
fig = figure()
ax = fig.add_subplot(111)
for fileId in files:
theta,dndtheta2=radial(fileId)
ax.plot(theta,dndtheta2,drawstyle='steps-mid',label=fileId)
ax.legend(loc="best",title="%3d - %3d GeV"%(Eliminf,Elimsup))
ax.set_yscale('log')
ax.grid(b=True,which='major')
ax.set_xlabel("$\\theta^2$ [deg$^2$]")
ax.set_ylabel("$dN/d\\theta^2$ [deg$^{-2}$]")
show()