diff --git a/EBL-spectrums.py b/EBL-spectrums.py index a1d8b4a..f6a2dc6 100755 --- a/EBL-spectrums.py +++ b/EBL-spectrums.py @@ -16,40 +16,40 @@ lamb,lambdaI = np.loadtxt("EBL_files/lambdaI_Dominguez.dat",unpack=True,usecols= hv = h*c/(lamb*1e-4) # erg density = 1e-6*(4*pi/c)*lambdaI/(hv**2) /(erg_to_GeV*1e9) *(1+z)**3 hv = hv *(erg_to_GeV*1e9) -ax.plot(hv,density,"--b",label="Dominguez and Al") +ax.plot(hv,density*hv**2,"--b",label="Dominguez and Al") # ==== Kneiske and Doll - "best fit" ==== hv,density = np.loadtxt("EBL_files/n_bestfit10.dat",unpack=True,usecols=[0,1]) -ax.plot(hv,density,"--r",label="Kneiske and Doll - 'best fit'") +ax.plot(hv,density*hv**2,"--r",label="Kneiske and Doll - 'best fit'") # ==== Kneiske and Doll - "lower limit" ==== hv,density = np.loadtxt("EBL_files/n_lowerlimit10.dat",unpack=True,usecols=[0,179]) #hv,density =np.loadtxt("EBL_files/n_lowerlimit10.dat",unpack=True,usecols=[0,1]) density=density*(1+z)**3 -ax.plot(hv,density,"--g",label="Kneiske and Doll - 'lower limit'") +ax.plot(hv,density*hv**2,"--g",label="Kneiske and Doll - 'lower limit'") # ==== Fransceschini ==== hv,density = np.loadtxt("EBL_files/n_Fra.dat",unpack=True,usecols=[0,11]) #hv,density = np.loadtxt("EBL_files/n_Fra.dat",unpack=True,usecols=[0,1]) density=density*(1+z)**3 -ax.plot(hv,density,"--c",label="Fraceschini") +ax.plot(hv,density*hv**2,"--c",label="Fraceschini") # ==== Finke ==== hv,density = np.loadtxt("EBL_files/n_Finke.dat",unpack=True,usecols=[0,201]) #hv,density = np.loadtxt("EBL_files/n_Finke.dat",unpack=True,usecols=[0,1]) density=density*(1+z)**3 -ax.plot(hv,density,"--m",label="Finke and Al") +ax.plot(hv,density*hv**2,"--m",label="Finke and Al") # ==== Gilmore ==== hv,density = np.loadtxt("EBL_files/n_Gil.dat",unpack=True,usecols=[0,14]) #hv,density = np.loadtxt("EBL_files/n_Gil.dat",unpack=True,usecols=[0,1]) -ax.plot(hv,density,"--y",label="Gilmore and Al") +ax.plot(hv,density*hv**2,"--y",label="Gilmore and Al") ax.set_xscale('log') ax.set_yscale('log') ax.grid(b=True,which='major') ax.legend(loc="best",title="z = %.0f"%z) ax.set_xlabel("energy [eV]") -ax.set_ylabel("$n_{EBL}$ [$cm^{-3} eV^{-1}$]") +ax.set_ylabel("$n_{EBL}$ [$eV.cm^{-3}$]") plt.show() diff --git a/Modules/Analytic.py b/Modules/Analytic.py index 03dbcc4..03c211a 100644 --- a/Modules/Analytic.py +++ b/Modules/Analytic.py @@ -4,6 +4,7 @@ from Read import ReadE_source, ReadD_source, ReadEGMF # Value for analytic expression Ee_default=1e3 #GeV +Egamma_default=1 #GeV B_default=1e-17 #Gauss lambda_B_default=0.3 #Mpc @@ -21,12 +22,24 @@ def Analytic_theta(E_ic,fileId): delta=Dic0/(2*RL0)*((Esource/2)**2 /Ee2 -1) return abs(arcsin(delta_to_theta*sin(delta))*degre) -def Analytic_delay(theta, E0_source, distSource): - delta_ic = arcsin(lambda_gg(E0_source)/distSource*theta/degre) - #c_delta_t = lambda_gg(E0_source)*(1-cos(delta_ic))-distSource*(1-cos(theta)) - #return abs(c_delta_t*Mpc/c) - coef = lambda_gg(E0_source)/(2*c/Mpc) - return coef*delta_ic#**2 +def Analytic_delay_vs_theta(theta,fileId): + E0_source = ReadE_source(fileId)*1e-3 #TeV + distSource = ReadD_source(fileId) + delta_ic = arcsin(distSource/lambda_gg(E0_source)*sin(theta/degre)) + c_delta_t = lambda_gg(E0_source)*(1-cos(delta_ic)) - distSource*(1-cos(theta/degre)) + return c_delta_t *Mpc/c # sec. + +def Analytic_delay_vs_Egamma(Egamma, fileId): + Esource=ReadE_source(fileId) #GeV + Dsource=ReadD_source(fileId) #Mpc + B=ReadEGMF(fileId) #Gauss + RL0=RL(Esource/2,B) + Dic0=Dic(Esource/2) + Ee2 = Ee(Egamma)**2 + delta_ic = Dic0/(2*RL0)*((Esource/2)**2 /Ee2 -1) + theta = arcsin(lambda_gg(Esource)/Dsource*sin(delta_ic/degre)) + c_delta_t = lambda_gg(Esource)*(1-cos(delta_ic)) - Dsource*(1-cos(theta/degre)) + return c_delta_t *Mpc/c # sec. # Compton accumulation def ECompton_threshold(Compton_threshold = 0.005): @@ -42,6 +55,9 @@ def lambdaIC(): def Eic(Ee=Ee_default): return 4*Ecmb*Ee**2/(3*me**2)*1e3 #GeV +def Ee(Egamma=Egamma_default): + return sqrt((3*me**2)/(4*Ecmb)*1e-3 *Egamma) #GeV + def tIC(): return lambdaIC()/(c*yr/Mpc) #yr @@ -66,6 +82,6 @@ def Ethreshold_ic(Ee=Ee_default,B=B_default): def Ethreshold_gg(): return (me)**2/Eebl *1e-3 #GeV -def lambda_gg(Egamma=1e3): # Egamma (GeV) - return 0.8 /(Egamma*1e-3) *1e3#Mpc (from Durrer and Neronov 2013) +def lambda_gg(Egamma=1): # Egamma (GeV) + return 800. /(Egamma)#Mpc (from Durrer and Neronov 2013) diff --git a/Modules/Analytic.pyc b/Modules/Analytic.pyc index 9e0d698..7fc8ea6 100644 Binary files a/Modules/Analytic.pyc and b/Modules/Analytic.pyc differ diff --git a/Modules/Angle.py b/Modules/Angle.py index f5fe33c..df0df32 100644 --- a/Modules/Angle.py +++ b/Modules/Angle.py @@ -44,11 +44,12 @@ def angle_vs_energy(fileName): nbBins = 40 energy =log10(energy) - angle,ener =histogram(energy,nbBins,weights=weight*theta) + angle,ener =histogram(energy,nbBins,weights=theta*weight) nb,ener =histogram(energy,nbBins,weights=weight) ener = 10**ener enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2 - angle = angle/nb + maxflux = ReadNphot_source(fileName) + angle = angle/(nb*maxflux) return enercenter, angle diff --git a/Modules/Angle.pyc b/Modules/Angle.pyc index 8e9cb79..0f5f0d9 100644 Binary files a/Modules/Angle.pyc and b/Modules/Angle.pyc differ diff --git a/Modules/Spectrum.pyc b/Modules/Spectrum.pyc index a95375e..fefe8f4 100644 Binary files a/Modules/Spectrum.pyc and b/Modules/Spectrum.pyc differ diff --git a/Modules/Timing.py b/Modules/Timing.py index 472205d..57dec87 100644 --- a/Modules/Timing.py +++ b/Modules/Timing.py @@ -2,10 +2,10 @@ from numpy import histogram, log10, shape, sqrt, arccos from matplotlib.pyplot import figure, show from scipy.integrate import quad from Read import ReadTime, ReadWeight, ReadGeneration, ReadE_source, ReadD_source -from Read import ReadNphot_source,ReadPosition,ReadMomentum,ReadNbLeptonsProd +from Read import ReadNphot_source,ReadPosition,ReadMomentum,ReadNbLeptonsProd,ReadEnergy from Constants import yr, degre from Integrand import comobileTime -from Analytic import Analytic_delay +from Analytic import Analytic_delay_vs_theta,Analytic_delay_vs_Egamma def timing(fileId,gen=-1): ''' @@ -93,8 +93,10 @@ def delay_vs_angle(fileId,gen=-1): ''' position = ReadPosition(fileId) momentum = ReadMomentum(fileId) + energy = ReadEnergy(fileId) weight = ReadWeight(fileId) time=ReadTime(fileId) + generation = ReadGeneration(fileId) hyp=sqrt(position[0]**2+position[1]**2+position[2]**2) position /= hyp @@ -103,21 +105,20 @@ def delay_vs_angle(fileId,gen=-1): costheta=position[0]*momentum[0]+position[1]*momentum[1]+position[2]*momentum[2] if gen != -1 : - generation = ReadGeneration(fileId) - cond=(abs(costheta)<=1) & (generation==gen) & (time>0) + cond=(abs(costheta)<=1) & (generation==gen) & (time>0) & (energy>1) else: - cond=(abs(costheta)<=1) & (time>0) + cond=(abs(costheta)<=1) & (time>0) & (energy>1) prop = float(shape(time[time<0])[0])/shape(time)[0] + print "gen", int(gen), "->", shape(time[time<0])[0], "negative time events among ",shape(time)[0], "~", prop theta = arccos(costheta[cond])*degre weight= weight[cond] time=time[cond]*yr - print "gen", int(gen), "->", shape(time)[0], "events", "negative time:",shape(time[time<0])[0], "~", prop cond=(theta!=0) - return theta[cond], time[cond], weight[cond] + return theta[cond], time[cond], weight[cond], generation[cond] def delay_vs_angle_histo(fileId,gen=-1): - theta, time, weight = delay_vs_angle(fileId,gen) + theta, time, weight, generation = delay_vs_angle(fileId,gen) nbBins = 100 theta=log10(theta) @@ -138,22 +139,17 @@ def drawDelay_vs_angle(fileId): fig = figure() ax = fig.add_subplot(111) - Nmax = ReadNphot_source(fileId) - - for gen in list(set(ReadGeneration(fileId))): - angle,delay,weight = delay_vs_angle(fileId,gen) - ax.plot(angle,delay/weight/Nmax,".",label="gen="+str(int(gen))) - - #angle,delay = delay_vs_angle_histo(fileId,gen) - #ax.plot(angle,delay,linestyle="steps-mid",label="gen="+str(int(gen))) + angle,delay,weight,generation = delay_vs_angle(fileId) + #for gen in list(set(generation)): + for gen in [2]: + cond= (generation==gen) + ax.plot(angle[cond],delay[cond]/weight[cond],".",label="gen="+str(int(gen))) angle,delay = delay_vs_angle_histo(fileId) - ax.plot(angle,delay/Nmax,color='k',linestyle="steps-mid",label="gen=all") + #ax.plot(angle,delay,color='k',linestyle="steps-mid",label="gen=all") - Esource = ReadE_source(fileId) - dSource = ReadD_source(fileId) - yfit = Analytic_delay(angle,Esource,dSource) - ax.plot(angle,yfit,'--k',linewidth=2) + delay_fit = Analytic_delay_vs_theta(angle,fileId) + ax.plot(angle,delay_fit,'--r',linewidth=2) ax.set_xscale('log') ax.set_yscale('log') @@ -163,3 +159,68 @@ def drawDelay_vs_angle(fileId): ax.set_ylabel("Time delay [s]") show() + +######## delay versus angle ############### +def delay_vs_energy(fileId,gen=-1): + ''' + Compute time delay versus arrival angle + Input: directory name + Input (optional): generation (default = all) + Output: arrival angle [degre], time delay [sec] + ''' + energy = ReadEnergy(fileId) + weight = ReadWeight(fileId) + time=ReadTime(fileId) + generation = ReadGeneration(fileId) + + if gen != -1 : + cond= (generation==gen) & (time>0) & (energy>1) + else: + cond= (time>0) & (energy>1) + + prop = float(shape(time[time<0])[0])/shape(time)[0] + print "gen", int(gen), "->", shape(time[time<0])[0], "negative time events among ",shape(time)[0], "~", prop + return energy[cond], time[cond]*yr, weight[cond], generation[cond] + +def delay_vs_energy_histo(fileId,gen=-1): + energy, time, weight, generation = delay_vs_energy(fileId,gen) + nbBins = 100 + + energy=log10(energy) + dt,dE=histogram(energy,nbBins,weights=time) + dN,dE=histogram(energy,nbBins,weights=weight) + dE=10**dE + Ecenter=(dE[1:nbBins+1]+dE[0:nbBins])/2 + timing = dt/dN + + return Ecenter, timing + +def drawDelay_vs_energy(fileId): + ''' + Plot angle versus time delay, generation by generation + Input: directory name + Output: graph of angle versus time delay, generation by generation + ''' + fig = figure() + ax = fig.add_subplot(111) + + energy,delay,weight,generation = delay_vs_energy(fileId) + #for gen in list(set(generation)): + for gen in [2]: + cond= (generation==gen) + ax.plot(energy[cond],delay[cond]/weight[cond],".",label="gen="+str(int(gen))) + + energy,delay = delay_vs_energy_histo(fileId) + #ax.plot(energy,delay,color='k',linestyle="steps-mid",label="gen=all") + + delay_fit = Analytic_delay_vs_Egamma(energy,fileId) + ax.plot(energy,delay_fit,'--r',linewidth=2) + + ax.set_xscale('log') + ax.set_yscale('log') + ax.grid(b=True,which='major') + ax.legend(loc="best") + ax.set_xlabel("E [GeV]") + ax.set_ylabel("Time delay [s]") + + show() diff --git a/Modules/Timing.pyc b/Modules/Timing.pyc index 50eeae8..ca9ca62 100644 Binary files a/Modules/Timing.pyc and b/Modules/Timing.pyc differ -- libgit2 0.21.2