Add 2d-histogram to arrival direction

- count vs theta and phi

Add 2d-histogram to arrival direction
- count vs theta and phi
Thomas Fitoussi
1 parent ab1642d3
Showing 3 changed files with 46 additions and 48 deletions Show diff stats
Arrival_direction.py
Constantes.py
Constantes.pyc
@@ -3,6 +3,7 @@
 from sys import argv
 import numpy as np
 import matplotlib.pyplot as plt
+from matplotlib.colors import LogNorm
 from mpl_toolkits.mplot3d.axes3d import Axes3D
 if len(argv)<2:
@@ -11,77 +12,73 @@ if len(argv)&lt;2:
 # figure: position of the arriving particles in space -> size circle = energy
 #=============================================================================
-nbBins = 50
+nbBins = 100
 degre= 180/np.pi
 borne=5 # degree
-Elim=1e-6 # GeV
+Elim=1e-1 # GeV
-fig = plt.figure()
+fig1 = plt.figure()
+fig2 = plt.figure()
 nbplots = len(argv)-1
 ind = 1
 for fileId in argv[1:]:
    energy,thetaDir,phiDir = np.loadtxt("Results_EGMF"+fileId+"/Results_momentum",unpack=True,usecols=[0,4,5])
    thetaPos,phiPos = np.loadtxt("Results_EGMF"+fileId+"/Results_position",unpack=True,usecols=[4,5])
+   charge,weight=np.loadtxt("Results_EGMF"+fileId+"/Results_extra",unpack=True,usecols=[0,2])
-   plotid = 100+nbplots*10+ind
-   ax = fig.add_subplot(plotid)
-
-   Energy = energy[energy>Elim]
-   print "Cut at", Elim, "GeV:", np.shape(Energy)
-   thetaPos = thetaPos[energy>Elim]*degre
-   phiPos = phiPos[energy>Elim]*degre
-   thetaDir = thetaDir[energy>Elim]*degre
-   phiDir = phiDir[energy>Elim]*degre
+   thetaPos = thetaPos*degre
+   phiPos = phiPos*degre
+   thetaDir = thetaDir*degre
+   phiDir = phiDir*degre
    theta = thetaDir - thetaPos
    phi = phiDir -phiPos
-   ## random position of the detector
-   #a=np.random.rand(2,1)
-   #w=a[0,0]*2.-1.
-   #u=np.sqrt(1-w*w)*np.cos(2*np.pi*a[1,0])
-   #v=np.sqrt(1.-w*w)*np.sin(2.*np.pi*a[1,0])
-   ## 2D -> use only u and v
-   #theta_detector = u*180 #degre
-   #phi_detector = v*90 #degre
+   cond=(energy>Elim) & (abs(theta)<borne) & (abs(phi)<borne)
+   Energy = energy[cond]
+   weight=weight[cond]
+   theta=theta[cond]
+   phi=phi[cond]
-   #FOV = 0.1 #degre
-
-   #energy = energy[energy>Elim 
-   #                  and np.absolute(theta-theta_detector)<FOV 
-   #                  and np.absolute(phi-phi_detector)<FOV]
-   #print "Cut at", Elim, "GeV:", np.shape(Energy)
-   #theta = theta_detector - theta[energy>Elim 
-   #               and np.absolute(theta-theta_detector)<FOV
-   #               and np.absolute(phi-phi_detector)<FOV]
-   #phi = phi_detector - phi[energy>Elim 
-   #            and np.absolute(theta-theta_detector)<FOV 
-   #            and np.absolute(phi-phi_detector)<FOV]
-   #Energy=np.append((Energy,energy),axis=1)
-   #Theta=np.append((Theta,theta),axis=1)
-   #Phi=np.append((Phi,phi),axis=1)
+   print "Cut at", Elim, "GeV:", np.shape(Energy)
    Emax = max(Energy)
    Emin = min(Energy)
    area = np.pi * (15 * Energy/Emax)**2 +15
    colors = Energy/Emax
-   sax = ax.scatter(theta, phi, s=area, c=colors, alpha=0.3, marker='o')
+
+   plotid = 100+nbplots*10+ind
+   #ax = fig1.add_subplot(plotid)
+   #sax = ax.scatter(theta, phi, s=area, c=colors, alpha=0.3, marker='o')
+   #ax.set_xlim((-borne,borne))
+   #ax.set_ylim((-borne,borne))
+   #ax.set_xlabel("$\\theta$ in degree")
+   #ax.set_ylabel("$\\phi$ in degree")
+   #ax.set_title("EGMF: $10^{-%.0f"%float(fileId)+"}$Gauss")
+ 
+   #if Emin < 1e-3:
+   #   liminf = "%.0f"%(Emin*1e6)+"keV"
+   #elif Emin < 1:
+   #   liminf = "%.0f"%(Emin*1e3)+"MeV"
+   #else: 
+   #   liminf = "%.0f"%(Emin)+"GeV"
+   #limsup = "%.0f"%(Emax*1e-3)+"TeV"
+   #cbar = fig1.colorbar(sax, ticks=[0.01,1])
+   #cbar.ax.set_yticklabels([liminf,limsup])
+   
+   ###########################################
+
+   ax = fig2.add_subplot(plotid)
+   H,xedges,yedges,im = ax.hist2d(theta,phi,nbBins,weights=weight,norm=LogNorm())
+   extent = [yedges[0], yedges[-1], xedges[0], xedges[-1]]
+   cbar=fig2.colorbar(im, ax=ax)
+   cbar.ax.set_ylabel("counts")
    ax.set_xlim((-borne,borne))
    ax.set_ylim((-borne,borne))
    ax.set_xlabel("$\\theta$ in degree")
    ax.set_ylabel("$\\phi$ in degree")
    ax.set_title("EGMF: $10^{-%.0f"%float(fileId)+"}$Gauss")
-
-   if Emin < 1e-3:
-      liminf = "%.0f"%(Emin*1e6)+"keV"
-   elif Emin < 1:
-      liminf = "%.0f"%(Emin*1e3)+"MeV"
-   else: 
-      liminf = "%.0f"%(Emin)+"GeV"
-   limsup = "%.0f"%(Emax*1e-3)+"TeV"
-   cbar = fig.colorbar(sax, ticks=[0.01,1])
-   cbar.ax.set_yticklabels([liminf,limsup])
-
+ 
    ind=ind+1
 plt.show()
@@ -15,7 +15,7 @@ lambdaC=hb/(m*c)
 Mpc=(3.0856776e+16)*1e8 # Mpc to cm
 erg_to_GeV=1/(1.602176565e-3)
 degre=180/pi
-s_to_yr=3600*24*365.25
+yr=3600*24*365.25 # s
 # Cosmology
 a0=1.
@@ -51,6 +51,7 @@ lambda_B=0.3 #Mpc
 Dic= 3*(me*1e-6)**2/(4*sigmaT*rhoCMB*1e-9*Ee) /Mpc #Mpc
 lambdaIC=1/(nCMB*sigmaT*Mpc) #Mpc
 Eic= 4*Ecmb*Ee**2/(3*me**2)*1e3 #GeV 
+tIC=lambdaIC/(c*yr/Mpc) #yr
 # Larmor radius
 RL=(Ee/erg_to_GeV)/(e*B) /Mpc #Mpc