Correction on the theta^2 distribution and on the analytic expression of theta vs E

Thomas Fitoussi
1 parent 138898c8
Showing 5 changed files with 48 additions and 33 deletions Show diff stats
Analysis.py
Modules/Analytic.py
Modules/Analytic.pyc
Modules/Angle.py
Modules/Angle.pyc
@@ -6,7 +6,7 @@ from Modules.Read import ReadEnergy, ReadTime, ReadExtraFile, ReadProfile
 from Modules.Read import ReadMomentumAngle, ReadPositionAngle, resultDirectory
 from Modules.Spectrum import spectrum
 from Modules.Map import computeMap
-from Modules.Angle import angle_vs_energy
+from Modules.Angle import angle_vs_energy, radial
 from Modules.Timing import timing
 from Modules.Constants import degre
@@ -29,8 +29,9 @@ for fileId in argv[1:]:
    energy = ReadEnergy(fileId)
    weightini, generation, theta_arrival, Esource, dir_source = ReadExtraFile(fileId,[2,3,4,5,6])
    #weightini, generation, theta_arrival, Esource = ReadExtraFile(fileId,[2,3,4,5])
-   nbPhotonsEmitted=ReadProfile(fileId,[3]) 
-   weightini /= nbPhotonsEmitted
+   n_phot, n_lept = ReadProfile(fileId,[3,5]) 
+   ratio = n_phot#/n_lept
+   weightini /= ratio
    theta_arrival*= degre
    thetaDir,phiDir = ReadMomentumAngle(fileId)*degre
@@ -68,7 +69,7 @@ for fileId in argv[1:]:
       nbBins = 100
       # draw source spectrum
       Es=array(list(set(Esource)))
-      Ws= (Es/min(Es))**(1-powerlaw_index) /nbPhotonsEmitted
+      Ws= (Es/min(Es))**(1-powerlaw_index)  / ratio
       Es,Fs = spectrum(Es,Ws,nbBins=nbBins)
       Es=Es[:,newaxis]
       Fs=Fs[:,newaxis]
@@ -92,8 +93,10 @@ for fileId in argv[1:]:
          Spectrum = append(Spectrum,flux_0,axis=1)
       #  IMAGING, RADIAL DISTRIBUTION AND ANGLE VERSUS ENERGY ======================#
-      nbBins = 60
-      theta2,dndtheta2,ener,angle = angle_vs_energy(theta_arrival,energy,weight,nbBins)
+      nbBins = 50
+      cond = energy>1e-3
+      ener,angle = angle_vs_energy(theta_arrival[cond],energy[cond],weight[cond],nbBins)
+      theta2,dndtheta2 = radial(theta[cond],weight[cond],nbBins)
       theta2=theta2[:,newaxis]
       dndtheta2=dndtheta2[:,newaxis]
       ener=ener[:,newaxis]
@@ -128,7 +131,7 @@ for fileId in argv[1:]:
    Jet_Opening=[180,60,30] # degre (180 <=> isotrop)
    powerlaw_index = 2
    Elim = 1e-1 # GeV
-   thetalim = 50 # degre
+   thetalim = 20 # degre
    # apply source spectrum
    weight_source = (Esource/min(Esource))**(1-powerlaw_index)
@@ -160,7 +163,7 @@ for fileId in argv[1:]:
       nbBins = 100
       # draw source spectrum
       Es=array(list(set(Esource)))
-      Ws= (Es/min(Es))**(1-powerlaw_index) /nbPhotonsEmitted
+      Ws= (Es/min(Es))**(1-powerlaw_index) / ratio
       Es,Fs = spectrum(Es,Ws,nbBins=nbBins)
       Es=Es[:,newaxis]
       Fs=Fs[:,newaxis]
@@ -184,9 +187,10 @@ for fileId in argv[1:]:
          Spectrum = append(Spectrum,flux_0,axis=1)
       #  IMAGING, RADIAL DISTRIBUTION AND ANGLE VERSUS ENERGY ======================#
-      nbBins = 60
+      nbBins = 50
       computeMap(theta,phi,weight,energy,fileId,jet_opening_angle,nbBins,borne=[thetalim,thetalim])
-      theta2,dndtheta2,ener,angle = angle_vs_energy(theta_arrival,energy,weight,nbBins)
+      ener,angle = angle_vs_energy(theta_arrival,energy,weight,nbBins)
+      theta2,dndtheta2 = radial(theta_arrival,weight,nbBins)
       theta2=theta2[:,newaxis]
       dndtheta2=dndtheta2[:,newaxis]
       ener=ener[:,newaxis]
@@ -201,9 +205,7 @@ for fileId in argv[1:]:
          Radial = append(Radial,dndtheta2,axis=1)
          Angle_Energy = append(Angle_Energy,angle,axis=1)
-      #=============================================================================#
-      #  TIME DISTRIBUTION AND TIME DELAY VERSUS ANGLE
-      #=============================================================================#
+      #  TIME DISTRIBUTION AND TIME DELAY VERSUS ANGLE =============================#
       nbBins = 100
       delta_t,dNdt = timing(time,weight,nbBins)
       delta_t=delta_t[:,newaxis]
@@ -16,9 +16,8 @@ def Analytic_theta(E_ic,fileId):
    delta_to_theta=lambda_gg(Esource)/Dsource
    RL0=RL(Esource/2,B)
    Dic0=Dic(Esource/2)
-   Ee2= E_ic /Eic(1) #GeV^2 
-   delta=Dic0/(2*RL0)*((Esource/2)**2 /Ee2 -1)
-   return abs(arcsin(delta_to_theta*sin(delta)))
+   delta=Dic0/(2*RL0)*((Esource/(2*Ee(E_ic)))**2 -1)
+   return abs(arcsin(delta_to_theta*sin(delta)))*degre
 def Analytic_delay_vs_theta(theta,fileId):
    Esource,distSource=ReadProfile(fileId,[0,2]) #GeV #Mpc 
@@ -60,7 +59,7 @@ 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 
+   return me*sqrt((3*Egamma*1e-3 )/(4*Ecmb)) #GeV 
 def tIC():
    return lambdaIC()/(c*yr/Mpc) #yr
@@ -87,4 +86,4 @@ def Ethreshold_gg():
    return (me)**2/Eebl *1e-3 #GeV
 def lambda_gg(Egamma=1): # Egamma (GeV)
-   return 800. /(Egamma)#Mpc (from Durrer and Neronov 2013)
+   return 800. /(Egamma) *1e3#Mpc (from Durrer and Neronov 2013)
@@ -19,20 +19,22 @@ def angle_vs_energy(theta,energy,weight,nbBins=50):
    enercenter = (ener[1:nbBins+1]+ener[0:nbBins])/2
    angle /= nb
-   #  RADIAL DISTRIBUTION ========================================================#
+   return enercenter, angle
+
+def radial(theta,weight,nbBins=50):
    cond = theta!=0
-   #theta = 2*log10(theta[cond])
-   theta = theta[cond]**2
+   theta = 2*log10(theta[cond])
    weight = weight[cond]
-   #dn,angle2 = histogram(theta,nbBins,range=[2*log10(1e-6),2*log10(90)],weights=weight)
-   dn,angle2 = histogram(theta,nbBins,range=[0,20],weights=weight)
-   #angle2=10**angle2
+   dn,angle2 = histogram(theta,nbBins,range=[2*log10(1e-3),2*log10(25)],weights=weight)
+   angle2=10**angle2
+   #theta = theta**2
+   #dn,angle2 = histogram(theta,nbBins,range=[0,25],weights=weight)
    thetacenter = (angle2[1:nbBins+1]+angle2[0:nbBins])/2
    dtheta = angle2[1:nbBins+1]-angle2[0:nbBins]
    dndtheta = dn/dtheta
    dndtheta /= max(dndtheta)
-   return thetacenter, dndtheta, enercenter, angle
+   return thetacenter, dndtheta
 def drawRadial(files,all_source_spectrum=False,all_jets=False):
    '''
@@ -61,15 +63,15 @@ def drawRadial(files,all_source_spectrum=False,all_jets=False):
          ax.plot(theta,dndtheta2,drawstyle='steps-mid',label=fileId)
    ax.legend(loc="best")
-   #ax.set_xscale('log')
+   ax.set_xscale('log')
    ax.set_yscale('log')
    ax.grid(b=True,which='major')
    ax.set_xlabel("$\\theta^2$ [deg$^2$]")
-   ax.set_ylabel("$dN/d\\theta$ [deg$^{-1}$]")
+   ax.set_ylabel("$dN/d\\theta^2$ [arbitrary]")
    show()
-def drawAngle_vs_energy(files):
+def drawAngle_vs_energy(files,all_source_spectrum=False,all_jets=False,PlotAnalytic=False):
    '''
       Plot arrival angle versus energy + analytic expression
       Input: list of directories
@@ -79,18 +81,30 @@ def drawAngle_vs_energy(files):
    ax1 = fig.add_subplot(111)
    for fileId in files:
-      ener,angle = ReadAngleVsEnergy(fileId,[0,1])
-      p=ax1.plot(ener,angle,drawstyle='steps-mid',label=fileId)
+      if all_source_spectrum:
+         i=1
+         for powerlaw_index in [1,1.5,2,2.5]: 
+            ener,angle = ReadAngleVsEnergy(fileId,[0,i])
+            p=ax1.plot(ener,angle,drawstyle='steps-mid',label="p=%.1f"%powerlaw_index)
+      elif all_jets:
+         i=5
+         for jet_opening in ["isotrop","60 deg","30 deg"]: 
+            ener,angle = ReadAngleVsEnergy(fileId,[0,i])
+            p=ax1.plot(ener,angle,drawstyle='steps-mid',label=jet_opening)
+      else:
+         ener,angle = ReadAngleVsEnergy(fileId,[0,1])
+         p=ax1.plot(ener,angle,drawstyle='steps-mid',label=fileId)
-      yfit = Analytic_theta(ener,fileId)
-      ax1.plot(ener,yfit,'--',color="m",linewidth=2,label="analytic")
+      if PlotAnalytic:
+         yfit = Analytic_theta(ener,fileId)
+         ax1.plot(ener,yfit,'--',color=p[0].get_color(),linewidth=2,label="analytic")
    ax1.legend(loc="best")
    ax1.set_xscale('log')
    ax1.set_yscale('log')
    ax1.set_ylim([0,180])
    ax1.grid(b=True,which='major')
-   ax1.set_ylabel("$\\theta$ [rad]")
+   ax1.set_ylabel("$\\theta$ [deg]")
    ax1.set_xlabel("energy [GeV]")
    show()