Time distribution

add a comparison between different precision on zlim Add a test on generation and charge in all scripts

Time distribution
add a comparison between different precision on zlim Add a test on generation and charge in all scripts
Thomas Fitoussi
1 parent d33c0beb
Showing 5 changed files with 55 additions and 30 deletions Show diff stats
Angle_distribution.py
Arrival_direction.py
Energy_distribution.py
Generation_distribution.py
Time_distribution.py
@@ -48,7 +48,7 @@ ind = 0
 for fileId in argv[2:]:
    energy,dirx,diry,dirz = loadtxt(fileName+fileId+"/Results_momentum",unpack=True,usecols=[0,1,2,3])
    posx,posy,posz = loadtxt(fileName+fileId+"/Results_position",unpack=True,usecols=[1,2,3])
-   weight,gen=loadtxt(fileName+fileId+"/Results_extra",unpack=True,usecols=[2,3])
+   charge,weight,gen=loadtxt(fileName+fileId+"/Results_extra",unpack=True,usecols=[0,2,3])
    hyp=sqrt(posx**2+posy**2+posz**2)
    posx=posx/hyp
@@ -59,7 +59,7 @@ for fileId in argv[2:]:
    diry=diry/hyp
    dirz=dirz/hyp
    costheta=dirx*posx+diry*posy+dirz*posz
-   cond=(costheta<=1)&(costheta>=-1)&(gen==2)
+   cond=(costheta<=1)&(costheta>=-1)&(charge==0)
    theta = arccos(costheta[cond])*degre
    weight= weight[cond]
    energy= energy[cond]
@@ -14,7 +14,7 @@ if len(argv)&lt;2:
 nbBins = 50
 degre= 180/np.pi
 borne=5 # degree
-Elim=1e-3 # GeV
+Elim=1e-6 # GeV
 fig = plt.figure()
@@ -39,7 +39,11 @@ Emin = min(Edata)
 ind = 0
 for fileId in argv[2:]:
    energy= np.loadtxt(fileName+fileId+"/Results_momentum",unpack=True,usecols=[0])
-   weight= np.loadtxt(fileName+fileId+"/Results_extra",unpack=True,usecols=[2])
+   charge,weight=np.loadtxt(fileName+fileId+"/Results_extra",unpack=True,usecols=[0,2])
+
+   cond= (charge==0)
+   energy=energy[cond]
+   weight=weight[cond]
    print "data shape:", np.shape(energy) 
    Edata =np.log10(energy)
@@ -54,14 +58,13 @@ for fileId in argv[2:]:
             label="$J_\gamma$ $10^{-%.0f"%float(fileId)+"}$Gauss")
    elif argv[1] == "EBL":
       ax.plot(enercenter,flux[:]/max(flux),"."+color[ind],label=Model[int(fileId)-1])
-
-   # Results from Elmag
-   Elmag=fileName+fileId+"/spec_diff_test"
-   if os.path.isfile(Elmag):
-      energy,fluxGamma,fluxLepton = np.loadtxt(Elmag,unpack=True,usecols=[0,1,2])
-      energy=energy*10**(-9)
-      flux=fluxGamma+fluxLepton
-      ax.plot(energy,flux/max(flux),"-"+color[ind],label="Elmag - "+Model[int(fileId)-1])
+      # Results from Elmag
+      Elmag=fileName+fileId+"/spec_diff_test"
+      if os.path.isfile(Elmag):
+         energy,fluxGamma,fluxLepton = np.loadtxt(Elmag,unpack=True,usecols=[0,1,2])
+         energy=energy*10**(-9)
+         flux=fluxGamma+fluxLepton
+         ax.plot(energy,flux/max(flux),"-"+color[ind],label="Elmag - "+Model[int(fileId)-1])
    ind=ind+1
@@ -10,7 +10,7 @@ if len(argv)&lt;2:
 # figure: energy distribution
 #=============================
-nbBins = 8 
+nbBins = 7
 color=['b','g','r','c','m','y']
 fig = plt.figure()
@@ -20,7 +20,7 @@ ind = 0
 for fileId in argv[1:]:
    weight,nbGen=np.loadtxt("Results_EGMF"+fileId+"/Results_extra",unpack=True,usecols=[2,3])
-   plt.hist(nbGen,nbBins, weights=weight,range=[1,9],log=1,facecolor=color[ind],alpha=.5,
+   plt.hist(nbGen,nbBins,range=[2,9],log=1,facecolor=color[ind],alpha=.5, weights=weight,
          label="$10^{-%.0f"%float(fileId)+"}$Gauss")
    ind=ind+1
@@ -2,10 +2,9 @@
 from sys import argv
 import os.path
-import numpy as np
-from scipy.optimize import curve_fit
 import matplotlib.pyplot as plt
-from Constantes import s_to_yr
+from Constantes import *
+from scipy.integrate import quad
 if argv[1] == "EGMF":
    fileName="Results_EGMF"
@@ -14,6 +13,9 @@ elif argv[1] == &quot;EBL&quot;:
    B=10**(-15)
    Model=["Kneiske and Doll - 'best fit'","Kneiske and Doll - 'lower limit'",
          "Franceschini","Finke and Al","Gilmore and Al","Dominguez and Al"]
+elif argv[1] == "prec":
+   fileName="Results_prec"
+   B=10**(-17)
 else:
    print "Used: ./Energy_distribution.py arg1 ind1 ind2 ..."
    print "        arg1 = EGMF or EBL "
@@ -25,39 +27,57 @@ else:
 # figure: energy distribution
 #=============================
-nbBins = 200
-convert= 180/np.pi
+zSource = 0.0308
+nbBins = 100
+convert= 180/pi
 color=['b','r','g','c','m','y']
 fig = plt.figure()
 ax = fig.add_subplot(111)
-time= np.loadtxt(fileName+argv[2]+"/Results_position",unpack=True,usecols=[0])
+def properTime(z):
+   return -1/(H0*(1+z)*sqrt(omegaM*(1+z)**3+omegaK*(1+z)**2+omegaL))
+
+tlim = quad(properTime,zSource,0)[0]/s_to_yr
+print tlim
+
+time= loadtxt(fileName+argv[2]+"/Results_position",unpack=True,usecols=[0])
+time=time[time>0]
+time=time/tlim
 tmax = max(time)
 tmin = min(time) 
 print tmin, tmax
 ind = 0
 for fileId in argv[2:]:
-   time= np.loadtxt(fileName+fileId+"/Results_position",unpack=True,usecols=[0])
-   weight,gen=np.loadtxt(fileName+fileId+"/Results_extra",unpack=True,usecols=[2,3])
+   time= loadtxt(fileName+fileId+"/Results_position",unpack=True,usecols=[0])
+   charge,weight,gen=loadtxt(fileName+fileId+"/Results_extra",unpack=True,usecols=[0,2,3])
+   time=time/tlim
+
+   prop = float(shape(time[time<0])[0])/shape(time)[0]
+   print "data shape",fileName+fileId,":", shape(time), "negative time:", shape(time[time<0]), "~", prop
-   cond=(gen<8)
+   cond= (gen<10) & (time>0) & (charge==0)
    time=time[cond]
    weight=weight[cond]
+   time=log10(time)
-   print "data shape",fileName+fileId,":", np.shape(time) 
-
-   dN,dt=np.histogram(time,nbBins,range=(tmin,tmax),weights=weight)
+   dN,dt=histogram(time,nbBins,range=(log10(tmin),log10(tmax)),weights=weight)
+   #time=10**time
    timecenter=(dt[1:nbBins+1]+dt[0:nbBins])/2
    binSize=dt[1:nbBins+1]-dt[0:nbBins]
    dNdt=dN/binSize
    if argv[1] == "EGMF":
-      #ax.hist(time,nbBins,weights=weight,range=(tmin,tmax),log=1,facecolor=color[ind],alpha=.5,
-      #   label="$10^{-%.0f"%float(fileId)+"}$Gauss")
+      ax.hist(time,nbBins,weights=weight,range=(log10(tmin),log10(tmax)),log=1,facecolor=color[ind],
+            alpha=.5, label="$10^{-%.0f"%float(fileId)+"}$Gauss")
       ax.plot(timecenter[dNdt!=0],dNdt[dNdt!=0],"-"+color[ind],
          label="$10^{-%.0f"%float(fileId)+"}$Gauss")
+   elif argv[1] == "prec":
+      ax.hist(time,nbBins,weights=weight,range=(log10(tmin),log10(tmax)),log=1,facecolor=color[ind],
+            alpha=.5)#, label="prec = $10^{-%.0f"%float(fileId)+"}$")
+      ax.plot(timecenter[dNdt!=0],dNdt[dNdt!=0],"-"+color[ind],
+            label="prec = $10^{-%.0f"%float(fileId)+"}$")
    elif argv[1] == "EBL":
       #plt.hist(time,nbBins,weights=weight,range=(tmin,tmax),log=1,facecolor=color[ind],alpha=.5,
       ax.plot(timecenter,dNdt,"."+color[ind],
@@ -69,7 +89,9 @@ for fileId in argv[2:]:
 #ax.set_yscale('log')
 ax.grid(b=True,which='major')
 ax.legend(loc="best")
-ax.set_xlabel("Time [yr]")
-ax.set_ylabel("$dN/dt$ [yr$^{-1}$]")
+if argv[1]=="prec":
+   ax.legend(loc="best",title="$10^{-17}$Gauss")
+ax.set_xlabel("Time [$\Delta t / t$ log scale]")
+ax.set_ylabel("$t\ dN/dt$")
 plt.show()