Add usefull integrand in Constantes.py

Thomas Fitoussi
1 parent 4448ea24
Showing 5 changed files with 34 additions and 23 deletions Show diff stats
Angle_distribution.py
Constantes.py
Constantes.pyc
Energy_distribution.py
Time_distribution.py
@@ -33,8 +33,8 @@ def fit_theta(E,B):
    delta=Dic0/(2*RL0)*((Esource/2)**2 *Eic/E-1)
    return abs(arcsin(delta_to_theta*sin(delta))*degre)
-#fig2 = plt.figure()
-#ax2 = fig2.add_subplot(111)
+fig2 = plt.figure()
+ax2 = fig2.add_subplot(111)
 fig1 = plt.figure()
 gs = gridspec.GridSpec(2, 1, height_ratios=[4,1]) 
@@ -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)&(charge==0)
+   cond=(costheta<=1)&(costheta>=-1)
    theta = arccos(costheta[cond])*degre
    weight= weight[cond]
    energy= energy[cond]
@@ -80,26 +80,29 @@ for fileId in argv[2:]:
       #ax11.plot(energy,theta,'.'+color[ind],label="$10^{-%.0f"%float(fileId)+"}$Gauss - MC")
       B=10**(-float(fileId))
       yfit = fit_theta(enercenter,B)
-      ax11.plot(enercenter,yfit,'--'+color[ind+1],linewidth=2,label="$10^{-%.0f"%float(fileId)+"}$Gauss - analytic")
+      ax11.plot(enercenter,yfit,'--'+color[ind],linewidth=2,label="$10^{-%.0f"%float(fileId)+"}$Gauss - analytic")
    elif argv[1] == "EBL":
       ax11.plot(enercenter,angle,'.'+color[ind],label=Model[int(fileId)-1])
       if ind==0: 
          yfit = fit_theta(enercenter,B)
-         ax11.plot(enercenter,yfit,'--'+color[ind+1],linewidth=2,label="$10^{-15}$Gauss - analytic")
+         ax11.plot(enercenter,yfit,'--'+color[ind],linewidth=2,label="$10^{-15}$Gauss - analytic")
    error=angle/yfit-1
    ax12.plot(enercenter,error,'+'+color[ind])
    # figure: radial distribution
    #=============================
-   #nbBins = 20
-   #theta = log10(theta)
-   #nb,angle1=histogram(theta, nbBins, weights=weight)
+   nbBins = 20
+   theta = theta[theta!=0]
+   weight=weight[theta!=0]
+   #theta = log10(theta**2)
+   nb,angle1=histogram(theta**2,nbBins,  range=[0,1],weights=weight)
    #angle1=10**angle1
-   #anglecenter=(angle1[1:nbBins+1]+angle1[0:nbBins])/2
-   #binSize=angle1[1:nbBins+1]-angle1[0:nbBins]
-   #ax2.plot(anglecenter,nb,'-'+color[ind],label="$10^{-%.0f"%float(fileId)+"}$Gauss - MC")
-   #ax2.hist(theta, nbBins, weights=weight,log=True,
+   anglecenter=(angle1[1:nbBins+1]+angle1[0:nbBins])/2
+   binSize=angle1[1:nbBins+1]-angle1[0:nbBins]
+   ax2.plot(anglecenter,nb,'-'+color[ind],
+         label="$10^{-%.0f"%float(fileId)+"}$Gauss - MC")
+   #ax2.hist(theta, nbBins, range=[0,0.5], weights=weight,log=True,
    #      facecolor=color[ind],alpha=.5,label="$10^{-%.0f"%float(fileId)+"}$Gauss")
    ind=ind+1
@@ -119,11 +122,11 @@ plt.ylim(min(error)-step, max(error)+step)
 xticklabels = ax11.get_xticklabels()
 plt.setp(xticklabels, visible=False)
-#ax2.legend(loc="best")
+ax2.legend(loc="best")
 #ax2.set_xscale('log')
-#ax2.set_yscale('log')
-#ax2.grid(b=True,which='major')
-#ax2.set_xlabel("$\\theta$ [degre]")
-#ax2.set_ylabel("Nb events * weight")
+ax2.set_yscale('log')
+ax2.grid(b=True,which='major')
+ax2.set_xlabel("$\\theta^2$ [degre]")
+ax2.set_ylabel("Nb events * weight")
 plt.show()
@@ -64,3 +64,14 @@ Ethreshold_ic = Eic*Dic/RL # GeV
 # Pair production
 Ethreshold_gg=(me)**2/Eebl *1e-3 #GeV
 lambda_gg=0.8 # (E_gamma/1TeV)^-1 Gpc (from Durrer and Neronov 2013)
+
+# usefull integrand
+def comobileTime(z):
+   return -1/(H0*(1+z)*sqrt(omegaM*(1+z)**3+omegaK*(1+z)**2+omegaL))
+
+def distPhoton(z):
+   return -c/(H0*a0*sqrt(omegaM*(1+z)**3+omegaK*(1+z)**2+omegaL))
+
+def distLepton(z,E):
+   beta = sqrt(1-m**2*c**4/E**2)
+   return -beta*c/(H0*a0*sqrt(omegaM*(1+z)**3+omegaK*(1+z)**2+omegaL))
@@ -24,7 +24,7 @@ else:
 # figure: energy distribution
 #=============================
-nbBins = 70
+nbBins = 100
 convert= 180/np.pi
 color=['b','r','g','c','m','y']
@@ -70,7 +70,7 @@ for fileId in argv[2:]:
 def fit(x,a,b,c):
    return a* x**b +c
-a=1
+a=1/2.5
 b=0.5
 c=0
@@ -35,10 +35,7 @@ color=[&#39;b&#39;,&#39;r&#39;,&#39;g&#39;,&#39;c&#39;,&#39;m&#39;,&#39;y&#39;]
 fig = plt.figure()
 ax = fig.add_subplot(111)
-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
+tlim = quad(comobileTime,zSource,0)[0]/s_to_yr
 print tlim
 time= loadtxt(fileName+argv[2]+"/Results_position",unpack=True,usecols=[0])