Spectrum.py
2.6 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
#!/usr/bin/python
from sys import argv
import matplotlib.pyplot as plt
from Constantes import *
if argv[1] == "EGMF":
fileName="Results_EGMF"
elif argv[1] == "EBL":
fileName="Results_EBL"
B=10**(-15)
Model=["Kneiske and Doll - 'best fit'","Kneiske and Doll - 'lower limit'",
"Franceschini","Finke and Al","Gilmore and Al","Dominguez and Al"]
else:
print "Used: ./Energy_distribution.py arg1 ind1 ind2 ..."
print " arg1 = EGMF or EBL "
print " ind1 = file number \n"
print "Examples: "
print " to study EGMF14: ./Energy_distribution.py EGMF 14 "
print " to study EBL1 and EBL2: ./Energy_distribution.py EBL 1 2 \n"
exit()
# figure: energy distribution
#=============================
nbBins = 100
convert= 180/pi
color=['b','r','g','c','m','y']
fig = plt.figure()
ax = fig.add_subplot(111)
energy= loadtxt(fileName+argv[2]+"/Results_momentum",unpack=True,usecols=[0])
Edata =log10(energy)
Emax = max(Edata)
Emin = min(Edata)
ind = 0
for fileId in argv[2:]:
energy= loadtxt(fileName+fileId+"/Results_momentum",unpack=True,usecols=[0])
charge,weight=loadtxt(fileName+fileId+"/Results_extra",unpack=True,usecols=[0,2])
cond= (charge==0)
energy=energy[cond]
weight=weight[cond]
print "data shape:", shape(energy)
Edata =log10(energy)
flux,ener=histogram(Edata,nbBins,range=(Emin,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
maxflux = shape(energy)[0]
if argv[1] == "EGMF":
ax.plot(enercenter,flux[:]/maxflux,color=color[ind],drawstyle='steps-mid',
label="$J_\gamma$ $10^{-%.0f"%float(fileId)+"}$Gauss")
elif argv[1] == "EBL":
ax.plot(enercenter,flux[:]/maxflux,color=color[ind],drawstyle='steps-mid',
label=Model[int(fileId)-1])
# Results from Elmag
Elmag=fileName+fileId+"/spec_diff_test"
if os.path.isfile(Elmag):
energy,fluxGamma,fluxLepton = loadtxt(Elmag,unpack=True,usecols=[0,1,2])
energy=energy*10**(-9)
flux=fluxGamma+fluxLepton
flux=flux*2.5e-12
ax.plot(energy,flux,"-"+color[ind],label="Elmag - "+Model[int(fileId)-1])
ind=ind+1
def fit(E_ic):
return me*1e3/4*sqrt(3e9/Ecmb)*1e-9*sqrt(E_ic) #GeV
xfit = logspace(log10(10**(-3)),log10(2*10**4),200)
yfit = fit(xfit)
ax.plot(xfit,yfit,'--m',linewidth=2,label="$E^{0.5}$")
ax.set_xscale('log')
ax.set_yscale('log')
ax.grid(b=True,which='major')
ax.legend(loc="best")
ax.set_xlabel("energy [GeV]")
ax.set_ylabel("$E^2 \\frac{dN}{dE}$ [GeV.photon]")
plt.show()