Add time delay versus energy and versus arrival angle

Thomas Fitoussi
1 parent f4246390
Showing 8 changed files with 116 additions and 38 deletions Show diff stats
EBL-spectrums.py
Modules/Analytic.py
Modules/Analytic.pyc
Modules/Angle.py
Modules/Angle.pyc
Modules/Spectrum.pyc
Modules/Timing.py
Modules/Timing.pyc
@@ -16,40 +16,40 @@ lamb,lambdaI = np.loadtxt(&quot;EBL_files/lambdaI_Dominguez.dat&quot;,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()
@@ -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)
  
@@ -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
  
@@ -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()