Map.py
2.3 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
from numpy import shape, pi
from matplotlib.pyplot import figure, show
from matplotlib.colors import LogNorm
from Read import ReadEnergy, ReadWeight, ReadMomentumAngle, ReadPositionAngle
from Constants import degre
Elim=1e-1 # GeV
borne=5 # degree
nbBins = 100
def thetaphi(fileId):
energy = ReadEnergy(fileId)
weight = ReadWeight(fileId)
thetaDir,phiDir = ReadMomentumAngle(fileId)*degre
thetaPos,phiPos = ReadPositionAngle(fileId)*degre
theta = thetaDir - thetaPos
phi = phiDir -phiPos
cond=(energy>Elim) & (abs(theta)<borne) & (abs(phi)<borne)
energy = energy[cond]
weight=weight[cond]
theta=theta[cond]
phi=phi[cond]
print "file", fileId, "->", shape(energy)[0], "events"
return energy, theta, phi, weight
def drawMap(files):
print "Cut at", Elim, "GeV"
fig = figure()
nbplots = len(files)
ind=1
for fileId in files:
energy, theta, phi, weight = thetaphi(fileId)
ax = fig.add_subplot(100+nbplots*10+ind)
Emax = max(energy)
Emin = min(energy)
colors = energy/Emax
area = pi * (15 * colors)**2 +15
sax = ax.scatter(theta, phi, c=colors, s=area, alpha=0.3, marker='o')
if Emin < 1e-3:
liminf = "%.0f"%(Emin*1e6)+"keV"
elif Emin < 1:
liminf = "%.0f"%(Emin*1e3)+"MeV"
else:
liminf = "%.0f"%(Emin)+"GeV"
limsup = "%.0f"%(Emax*1e-3)+"TeV"
cbar = fig.colorbar(sax, ticks=[0.01,1])
cbar.ax.set_yticklabels([liminf,limsup])
ax.set_xlim((-borne,borne))
ax.set_ylim((-borne,borne))
ax.set_xlabel("$\\theta$ [deg]")
ax.set_ylabel("$\\phi$ [deg]")
ax.set_title(fileId)
ind=ind+1
show()
def drawMapHisto(files):
print "Cut at", Elim, "GeV"
fig = figure()
nbplots = len(files)
ind=1
for fileId in files:
energy, theta, phi, weight = thetaphi(fileId)
ax = fig.add_subplot(100+nbplots*10+ind)
H,xedges,yedges,im = ax.hist2d(theta,phi,nbBins,weights=weight/shape(weight),norm=LogNorm())
extent = [yedges[0], yedges[-1], xedges[0], xedges[-1]]
cbar=fig.colorbar(im, ax=ax)
cbar.ax.set_ylabel("counts")
ax.set_xlim((-borne,borne))
ax.set_ylim((-borne,borne))
ax.set_xlabel("$\\theta$ [deg]")
ax.set_ylabel("$\\phi$ [deg]")
ax.set_title(fileId)
ind=ind+1
show()