compute_EGB.py
2.5 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
#!/bin/python
from numpy import loadtxt, interp, zeros, size, savetxt, log, trapz, array, pi, nan_to_num
from src.read import select_events, ReadProfile
from src.distribution import distribution
from src.analytic import distance, Mpc, MeV
EBL_models=["Dominguez"]
redshifts=["0.001168203","0.001561193","0.002086388","0.002788261","0.003726248",
"0.004979780","0.006655007","0.008893791","0.011885710","0.015884140",
"0.021227660","0.028368760","0.037912180","0.050666060","0.067710430",
"0.090488610","0.120929500","0.161610900","0.215977800","0.288634000",
"0.385732100","0.515494700","0.688910100","0.920663500","1.230380000",
"1.644287000","2.197435000","2.936665000","3.924576000","5.244827000"]
Emin = 0.1074077 # GeV
nE= 50
SaveResults=zeros((size(redshifts)+1,nE+1))
SaveResults[0,0] = 0
SaveResults0=zeros((size(redshifts)+1,nE+1))
SaveResults0[0,0] = 0
# Ajello et Al matrix
# > Emissivity (ph/MeV/cm^2/s/sr/bin-z) as a function of redshift and Observed Energy (MeV)
# > First line are Energies in MeV
# > First element of each subsequent line is the REDSHIFT
mat = loadtxt("Results/EmissMatrix__PLE_Paper_10TeV.txt")
Eajello = mat[0,1:]*1e-3 # MeV --> GeV
Emiss = nan_to_num(mat[1:,1:])*1e3 # ph/MeV/cm^2/s/sr/bin-z --> ph/GeV/cm^2/s/sr/bin-z
for EBL in EBL_models:
print "EBL model:", EBL
i=1
for z in redshifts:
fileId="Simulations/z="+z+"/"+EBL
Emax = ReadProfile(fileId,cols=[1])*1e3 # GeV
print "\t > redshift", z, "Emax =", Emax, "GeV", log(Emax/Emin)
weight, energy, time, theta, phi, Esource, gen = select_events(fileId,Erange=[Emin,Emax])
weight *= Esource*log(Emax/Emin) # such as int_Emin^Emax dNs/dE dE = 1
# recompute weight according to Ajello et Al matrix
dNdE_ajello = interp(Esource,Eajello*(1+float(z)),Emiss[i-1,:])/(1+float(z))
weight *= dNdE_ajello#*16*pi**2*(distance(float(z))[1]*Mpc)**2
E,dNdE = distribution(energy,weight,nbBins=nE,x_range=[Emin,110e3])
SaveResults[0,1:] = E
SaveResults[i,1:] = dNdE
SaveResults[i,0] = float(z)
cond = (gen==0)
E,dNdE = distribution(energy[cond],weight[cond],nbBins=nE,x_range=[Emin,110e3])
SaveResults0[0,1:] = E
SaveResults0[i,1:] = dNdE
SaveResults0[i,0] = float(z)
i+=1
saving_file = "Results/EGB_"+EBL+"-gen=0.dat"
savetxt(saving_file,SaveResults0,fmt='%1.4e')
saving_file = "Results/EGB_"+EBL+".dat"
savetxt(saving_file,SaveResults,fmt='%1.4e')
print "\t >>> Results saved in ", saving_file