Spectrum.py
2.88 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
from numpy import histogram, log10, sqrt, array
from matplotlib.pyplot import figure, show, setp
from matplotlib import gridspec
from Read import ReadSpectrum, ReadSourceSpectrum
from Constants import degre
from Analytic import Ethreshold_gg
def spectrum(energy,weight,generation=[],nbBins=100):
'''
Compute flux versus energy
Input: events energy and weight (optional: number of bins)
Output: energy, flux (E^2 dN/dE)
'''
Emin=min(energy)
Emax=max(energy)
energy =log10(energy)
flux,ener=histogram(energy,nbBins,range=[log10(Emin),log10(Emax)],weights=weight)
ener = 10**ener
enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2
binSize=ener[1:nbBins+1]-ener[0:nbBins]
flux = enercenter**2 * flux/binSize
if generation == []:
return enercenter,flux
else:
cond=(generation==0)
flux_0,ener_0=histogram(energy[cond],nbBins,range=[log10(Emin),log10(Emax)],weights=weight[cond])
ener_0 = 10**ener_0
enercenter_0=(ener_0[1:nbBins+1]+ener_0[0:nbBins])/2
binSize=ener_0[1:nbBins+1]-ener_0[0:nbBins]
flux_0 = enercenter_0**2 * flux_0/binSize
return enercenter,flux,flux_0
def drawSpectrum(files,PlotAnalytic=False):
'''
Plot flux versus energy
Input: list of directories
Input (optional, default=False): plot analytic expression
Output: graph spectrum normalize to 1 source photon
'''
fig = figure()
gs = gridspec.GridSpec(2, 1, height_ratios=[3,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:
i=1#3
j=1#5
for powerlaw_index in [1,1.5,2,2.5]:
# read files
Es,Fs = ReadSourceSpectrum(fileId,[0,i])
p = ax1.plot(Es,Fs,linestyle=':')
# draw full spectrum
ener,flux,flux_0 = ReadSpectrum(fileId,[0,j,j+1])
ax1.plot(ener,flux,color=p[0].get_color(),drawstyle='steps-mid',label="p=%.1f"%powerlaw_index)
# primary gamma-rays contribution
ax1.plot(ener,flux_0,color=p[0].get_color(),linestyle='--',drawstyle='steps-mid')
contrib=flux_0/flux
ax2.plot(ener,contrib,drawstyle='steps-mid',color=p[0].get_color())
i+=1
j+=2
# add spectrum for Fermi/LAT and CTA (FoV?)
if PlotAnalytic:
minEner=min(ener)
alpha=flux[ener==minEner]/sqrt(minEner)
ax1.plot(ener,alpha*sqrt(ener),'--m',linewidth=2)
ax1.plot(ener,ener**(0)*max(flux)*0.95,'--m',linewidth=2)
ax1.axvline(x=Ethreshold_gg(), ymin=0., ymax = 1., color='r', linewidth=2)
ax1.legend(loc="best")
ax1.set_xscale('log')
ax1.set_yscale('log')
ax1.grid(b=True,which='major')
ax1.set_ylabel("$E^2 dN/dE$ [GeV]")
ax2.set_xscale('log')
ax2.grid(b=True,which='major')
ax2.set_xlabel("energy [GeV]")
xticklabels = ax1.get_xticklabels()
setp(xticklabels, visible=False)
show()