Timing.py
1.33 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
from numpy import histogram, log10, shape
from matplotlib.pyplot import figure, show
from scipy.integrate import quad
from Read import ReadTime, ReadWeight, Readz_source, ReadNphot_source
from Constants import yr
from Integrand import comobileTime
def timing(fileName):
time=ReadTime(fileName)
weight=ReadWeight(fileName)
prop = float(shape(time[time<0])[0])/shape(time)[0]
print "file", fileName, "->", shape(time)[0], "events", "negative time:", shape(time[time<0])[0], "~", prop
time=time[time>0]
weight=weight[time>0]
nbBins = 100
zSource = Readz_source(fileName)
tlim = quad(comobileTime,zSource,0)[0]/yr
time = log10(time/tlim)
dN,dt=histogram(time,nbBins,weights=weight)
timecenter=(dt[1:nbBins+1]+dt[0:nbBins])/2
binSize=dt[1:nbBins+1]-dt[0:nbBins]
Nmax=ReadNphot_source(fileName) # Nb photons emitted by the source
dNdt=dN/(Nmax*binSize)
maxflux = ReadNphot_source(fileName)
dNdt /= maxflux
return timecenter, dNdt
def drawTiming(files):
fig = figure()
ax = fig.add_subplot(111)
for fileId in files:
time,dNdt = timing(fileId)
ax.plot(time[dNdt!=0],dNdt[dNdt!=0],linestyle="steps-mid",label=fileId)
ax.set_yscale('log')
ax.grid(b=True,which='major')
ax.legend(loc="best")
ax.set_xlabel("Time [$\Delta t / t$ log scale]")
ax.set_ylabel("$t\ dN/dt$")
show()