Improve angular distribution (plot in log-log)

Thomas Fitoussi
1 parent e6f5876a
Showing 16 changed files with 269 additions and 49 deletions Show diff stats
Analysis.py
Modules/Analytic.py
Modules/Analytic.pyc
Modules/Angle.py
Modules/Angle.pyc
Modules/Constants.pyc
Modules/Generation.py
Modules/Generation.pyc
Modules/Map.py
Modules/Map.pyc
Modules/Spectrum.py
Modules/Spectrum.pyc
Modules/Timing.py
Modules/Timing.pyc
Modules/Weight.py
Modules/Weight.pyc
@@ -26,7 +26,7 @@ if shape(argv)[0] &lt; 3:
    error("not enough argument")
  
 elif argv[1] == "spectrum": 
-   drawSpectrum(argv[2:])
+   drawSpectrum(argv[2:],PlotAnalytic=True)
  
 elif argv[1] == "spectrumLept": 
    drawSpectrum_normaLepton(argv[2:])
 from Constants import *
-from numpy import sqrt, abs, sin, arcsin
+from numpy import sqrt, abs, cos, sin, arcsin
 from Read import ReadE_source, ReadD_source, ReadEGMF
  
 # Value for analytic expression
@@ -21,6 +21,13 @@ 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
+
 # Compton accumulation
 def ECompton_threshold(Compton_threshold = 0.005):
    return Compton_threshold/(4/3*Ecmb/me*1e-3) *me*1e-6 #GeV
@@ -3,12 +3,18 @@ from matplotlib.pyplot import figure, show, ylim, setp
 from matplotlib import gridspec
 from Read import ReadMomentum, ReadPosition, ReadEnergy, ReadWeight, ReadNphot_source
 from Analytic import Analytic_theta
+from Map import thetaphi
 from Constants import degre
  
 Eliminf = 1e0 #GeV
 Elimsup = 1e5 #GeV
  
 def compute_theta(fileName):
+   '''
+      Compute arrival angle
+      Input: directory name
+      Output: energy, weight, arrival angle
+   '''
    position = ReadPosition(fileName)
    momentum = ReadMomentum(fileName)
    energy = ReadEnergy(fileName)
@@ -29,12 +35,17 @@ def compute_theta(fileName):
    return energy, weight, theta
  
 def angle_vs_energy(fileName):
+   '''
+      Compute arrival angle versus energy
+      Input: directory name
+      Output: energy, arrival angle
+   '''
    energy,weight,theta = compute_theta(fileName)
    nbBins = 40
  
    energy =log10(energy)
-   angle,ener =histogram(energy,nbBins,weights=theta)
-   nb,ener =histogram(energy,nbBins)
+   angle,ener =histogram(energy,nbBins,weights=weight*theta)
+   nb,ener =histogram(energy,nbBins,weights=weight)
    ener = 10**ener
    enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2
    angle = angle/nb
@@ -42,6 +53,11 @@ def angle_vs_energy(fileName):
    return enercenter, angle
  
 def drawAngle(files):
+   '''
+      Plot arrival angle versus energy + analytic expression
+      Input: list of directories
+      Output: graph of arrival angle versus energy
+   '''
    fig = figure()
    gs = gridspec.GridSpec(2, 1, height_ratios=[4,1]) 
    fig.subplots_adjust(hspace=0.001)
@@ -75,31 +91,49 @@ def drawAngle(files):
  
    show()
  
-def radial(fileName):
-   nbBins = 20
-   energy,weight,theta = compute_theta(fileName)
-
-   dn,angle=histogram(theta**2,nbBins,range=[0,0.25],weights=weight)
-   thetacenter=(angle[1:nbBins+1]+angle[0:nbBins])/2
-   dtheta=angle[1:nbBins+1]-angle[0:nbBins]
-   dndtheta2=dn/dtheta
-   maxflux = ReadNphot_source(fileName)
-   dndtheta2 /= maxflux
-
-   return thetacenter, dndtheta2
-
-def drawRadial(files):
+##############################################################################################
+def radial(fileId,thetamax=90):
+   '''
+      Compute flux versus arrival angle
+      Input: directory name
+      Input (optionnal): maximal angle desired (by default full sky)
+      Output: arrival angle, flux (dn/dtheta)
+   '''
+   nbBins = 50
+   energy,weight,theta = compute_theta(fileId)
+
+   theta=log10(theta)
+   dn,angle = histogram(theta,nbBins,range=[log10(1e-6),log10(thetamax)],weights=weight)
+   angle=10**angle
+   thetacenter = (angle[1:nbBins+1]+angle[0:nbBins])/2
+   dtheta = angle[1:nbBins+1]-angle[0:nbBins]
+   dndtheta = dn/dtheta
+   maxflux = ReadNphot_source(fileId)
+   dndtheta /= maxflux
+
+   print "integral: ", sum(dtheta*dndtheta)
+
+   return thetacenter, dndtheta
+
+def drawRadial(files,thetamax=90):
+   '''
+      Plot flux versus arrival angle
+      Input: list of directories
+      Input (optionnal): maximal angle desired (by default full sky)
+      Output: graph of flux versus angle
+   '''
    fig = figure()
    ax = fig.add_subplot(111)
  
    for fileId in files:
-      theta,dndtheta2=radial(fileId)
+      theta,dndtheta2=radial(fileId,thetamax)
       ax.plot(theta,dndtheta2,drawstyle='steps-mid',label=fileId)
  
-   ax.legend(loc="best",title="%3d - %3d GeV"%(Eliminf,Elimsup))
+   ax.legend(loc="best")
+   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^2$ [deg$^{-2}$]")
+   ax.set_xlabel("$\\theta$ [deg]")
+   ax.set_ylabel("$dN/d\\theta$ [deg$^{-1}$]")
  
    show()
@@ -4,6 +4,11 @@ from Read import ReadGeneration, ReadWeight
  
  
 def drawGeneration(files):
+   '''
+      Plot histogram of particles'generations
+      Input: list of directories
+      Output: histogram of generations
+   '''
    fig = figure()
    ax = fig.add_subplot(111)
    nbBins = 7
-from numpy import shape, pi, histogram2d, arange, newaxis, exp, log, array
+from numpy import shape, pi, histogram2d, arange, newaxis, exp, log, array, sin, cos
 from scipy.signal import convolve2d
 from matplotlib.pyplot import figure, show, matshow
 from matplotlib.colors import LogNorm
@@ -7,14 +7,21 @@ from Read import ReadD_source, Readz_source, ReadEGMF
 from Constants import degre
  
 thetalim = 180
-philim = 90
-step = 1 #deg
+philim = 180
+step = 0.25 #deg
 nbBins = array([int(2*philim/step)+1,int(2*thetalim/step)+1])
 philim=float(philim)
 thetalim=float(thetalim)
 limits = [[-thetalim,thetalim],[-philim,philim]]
  
 def thetaphi(fileId,Elim=0.1):
+   '''
+      Compute arrival angles (theta,phi) selected with energy upper than Elim 
+      normalize to 1 source photon
+      Input: directory name
+      Input (optional): lower energy limit in GeV (default = 100MeV)
+      Output: energy (GeV), theta, phi (degre), weight
+   '''
    energy = ReadEnergy(fileId)
    weight = ReadWeight(fileId)
    thetaDir,phiDir = ReadMomentumAngle(fileId)*degre
@@ -42,7 +49,15 @@ def Gaussian2D(center=[0,0],fwhm=0.01):
    y = y[:,newaxis]
    return exp(-4*log(2) * ((x-center[0])**2 + (y+center[1])**2) / fwhm**2)
  
-def drawMap(files,Elim=0.1,source=isotrop(),borne=[180,90]):
+def drawMap(files,Elim=0.1,source=isotrop(),borne=[180,180]):
+   '''
+      Plot 2D histogram of arrival direction of the photons with energy upper than Elim
+      Input: list of directories
+      Input (optional): lower energy limit (in GeV (default = 100MeV)
+      Input (optional): nature of the source (default isotrop) 
+      Input (optional): limits on theta, phi 
+      Output: 2D histogram of theta cos(phi), theta sin(phi)
+   '''
    print "Cut at", Elim, "GeV" 
    fig = figure()
    nbplots = len(files)
@@ -51,7 +66,7 @@ def drawMap(files,Elim=0.1,source=isotrop(),borne=[180,90]):
    for fileId in files:
       energy, theta, phi, weight = thetaphi(fileId,Elim)
  
-      H,xedges,yedges = histogram2d(phi,theta,nbBins,weights=weight,range=limits)
+      H,xedges,yedges = histogram2d(theta*cos(phi),theta*sin(phi),nbBins,weights=weight,range=limits)
       conv=convolve2d(H,source,boundary="wrap",mode="same")
  
       ax = fig.add_subplot(100+nbplots*10+ind)
@@ -62,8 +77,8 @@ def drawMap(files,Elim=0.1,source=isotrop(),borne=[180,90]):
       cbar.ax.set_ylabel("counts normalize to 1 photon emitted")
       ax.set_xlim((-borne[0],borne[0]))
       ax.set_ylim((-borne[1],borne[1]))
-      ax.set_xlabel("$\\theta$ [deg]")
-      ax.set_ylabel("$\\phi$ [deg]")
+      ax.set_xlabel("$\\theta$ $\\cos(\\phi)$ [deg]")
+      ax.set_ylabel("$\\theta$ $\\sin(\\phi)$ [deg]")
       ax.grid(b=True,which='major')
       zsource=Readz_source(fileId)
       Dsource=ReadD_source(fileId)
@@ -7,6 +7,11 @@ from Read import ReadNphot_source, ReadGeneration
 from Analytic import Analytic_flux, Ethreshold_gg
  
 def spectrum(fileId,gen=-1):
+   '''
+      Compute flux versus energy
+      Input: directory name, generation (optional, by default: all generation)
+      Output: energy, flux
+   '''
    energy = ReadEnergy(fileId)
    weight = ReadWeight(fileId)
  
@@ -30,12 +35,19 @@ def spectrum(fileId,gen=-1):
    return enercenter,flux
  
 ###########################################################################################
-def drawSpectrum(files,PlotElmag=False,PlotAnalytic=False):
+def drawSpectrum(files,gen=-1,PlotElmag=False,PlotAnalytic=False):
+   '''
+      Plot flux versus energy
+      Input: list of directories
+      Input (optional, default=all): selection on the generation
+      Input (optional, default=False): Plot Elmag spectrum if exists, plot analytic expression
+      Output: graph spectrum normalize to 1 source photon
+   '''
    fig = figure()
    ax = fig.add_subplot(111)
  
    for fileId in files:
-      energy,flux = spectrum(fileId)
+      energy,flux = spectrum(fileId,gen)
       if PlotElmag:
          maxflux = max(flux)*(.8)
       p=ax.plot(energy,flux,drawstyle='steps-mid',label=fileId)
@@ -47,9 +59,10 @@ def drawSpectrum(files,PlotElmag=False,PlotAnalytic=False):
          ax.plot(energy,flux,"--",color=p[0].get_color(),label="Elmag - "+fileId)
  
    if PlotAnalytic:
-      alpha=2.8e3
+      minEner=min(energy)
+      alpha=flux[energy==minEner]/sqrt(minEner)
       ax.plot(energy,alpha*sqrt(energy),'--m',linewidth=2)
-      ax.plot(energy,energy**(0)*max(flux),'--m',linewidth=2)
+      ax.plot(energy,energy**(0)*max(flux)*0.95,'--m',linewidth=2)
       ax.axvline(x=Ethreshold_gg(), ymin=0., ymax = 1., color='r', linewidth=2)
  
    ax.set_xscale('log')
@@ -65,6 +78,11 @@ def drawSpectrum(files,PlotElmag=False,PlotAnalytic=False):
    show()
  
 def drawSpectrum_normaLepton(files):
+   '''
+      Plot flux versus energy compare with analytic expression
+      Input: list of directories
+      Output: graph spectrum normalize to 1 lepton generated
+   '''
    fig = figure()
    gs = gridspec.GridSpec(2, 1, height_ratios=[4,1]) 
    fig.subplots_adjust(hspace=0.001)
@@ -101,6 +119,11 @@ def drawSpectrum_normaLepton(files):
  
 ########################################################################################
 def drawSpectrumGen(fileId):
+   '''
+      Plot flux versus energy generation by generation
+      Input: directory name
+      Output: graph spectrum normalize to 1 photon source
+   '''
    fig = figure()
    ax = fig.add_subplot(111)
  
@@ -110,8 +133,9 @@ def drawSpectrumGen(fileId):
  
    energy,flux = spectrum(fileId)
    p=ax.plot(energy,flux,drawstyle='steps-mid',color='k',label="gen = all")
-   alpha=2e3
  
+   minEner=min(energy)
+   alpha=flux[energy==minEner]/sqrt(minEner)
    ax.plot(energy,sqrt(energy)*alpha,'--m',linewidth=2)
    ax.plot(energy,energy**(0)*max(flux),'--m',linewidth=2)
    ax.axvline(x=Ethreshold_gg(), ymin=0., ymax = 1., color='r', linewidth=2)
-from numpy import histogram, log10, shape
+from numpy import histogram, log10, shape, sqrt, arccos
 from matplotlib.pyplot import figure, show
 from scipy.integrate import quad
-from Read import ReadTime, ReadWeight, Readz_source, ReadNphot_source
-from Constants import yr
+from Read import ReadTime, ReadWeight, ReadGeneration, ReadE_source, ReadD_source
+from Read import ReadNphot_source,ReadPosition,ReadMomentum,ReadNbLeptonsProd
+from Constants import yr, degre
 from Integrand import comobileTime
+from Analytic import Analytic_delay
  
-def timing(fileName):
-   time=ReadTime(fileName)
-   weight=ReadWeight(fileName)
-   prop = float(shape(time[time<0])[0])/shape(time)[0]
-   print "file", fileName, "->", shape(time)[0], "events", "negative time:", shape(time[time<0])[0], "~", prop
+def timing(fileId,gen=-1):
+   '''
+      Compute flux versus time delay
+      Input: directory name
+      Input (optional): generation (default = all)
+      Output: time delay (sec), flux
+   '''
+   time=ReadTime(fileId)
+   weight=ReadWeight(fileId)
+
+   if gen != -1:
+      generation = ReadGeneration(fileId)
+      time=time[generation==gen]
+      weight=weight[generation==gen]
+      prop = float(shape(time[time<0])[0])/shape(time)[0]
+      print "gen", int(gen), "->", shape(time)[0], "events", "negative time:",shape(time[time<0])[0], "~", prop
+   else:
+      prop = float(shape(time[time<0])[0])/shape(time)[0]
+      print "file", fileId, "->", shape(time)[0], "events", "negative time:",shape(time[time<0])[0], "~", prop
  
    time=time[time>0]
    weight=weight[time>0]
  
    nbBins = 100
-   zSource = Readz_source(fileName)
-   tlim = quad(comobileTime,zSource,0)[0]/yr
-   time = log10(time/tlim)
+   time = log10(time*yr)
  
    dN,dt=histogram(time,nbBins,weights=weight)
    timecenter=(dt[1:nbBins+1]+dt[0:nbBins])/2
    binSize=dt[1:nbBins+1]-dt[0:nbBins]
-   Nmax=ReadNphot_source(fileName) # Nb photons emitted by the source
+   Nmax=ReadNphot_source(fileId) # Nb photons emitted by the source
    dNdt=dN/(Nmax*binSize)
-   maxflux = ReadNphot_source(fileName)
+   maxflux = ReadNphot_source(fileId)
    dNdt /= maxflux
  
-   return timecenter, dNdt
+   return 10**timecenter, dNdt
  
 def drawTiming(files):
+   '''
+      Plot flux versus time delay
+      Input: list of directories
+      Output: graph of flux versus time delay 
+   '''
    fig = figure()
    ax = fig.add_subplot(111)
    for fileId in files:
       time,dNdt = timing(fileId) 
       ax.plot(time[dNdt!=0],dNdt[dNdt!=0],linestyle="steps-mid",label=fileId)
  
+   ax.set_xscale('log')
+   ax.set_yscale('log')
+   ax.grid(b=True,which='major')
+   ax.legend(loc="best")
+   ax.set_xlabel("Time delay [s]")
+   ax.set_ylabel("$t\ dN/dt$")
+
+   show()
+
+def drawTimingGen(fileId):
+   fig = figure()
+   ax = fig.add_subplot(111)
+   for gen in list(set(ReadGeneration(fileId))):
+      time,dNdt = timing(fileId,gen) 
+      ax.plot(time[dNdt!=0],dNdt[dNdt!=0],linestyle="steps-mid",label="gen="+str(int(gen)))
+
+   time,dNdt = timing(fileId) 
+   ax.plot(time[dNdt!=0],dNdt[dNdt!=0],color='k',linestyle="steps-mid",label="gen=all")
+
+   ax.set_xscale('log')
    ax.set_yscale('log')
    ax.grid(b=True,which='major')
    ax.legend(loc="best")
-   ax.set_xlabel("Time [$\Delta t / t$ log scale]")
+   ax.set_xlabel("Time delay [s]")
    ax.set_ylabel("$t\ dN/dt$")
  
    show()
+
+######## delay versus angle ###############
+def delay_vs_angle(fileId,gen=-1):
+   '''
+      Compute time delay versus arrival angle
+      Input: directory name
+      Input (optional): generation (default = all)
+      Output: arrival angle [degre], time delay [sec]
+   '''
+   position = ReadPosition(fileId)
+   momentum = ReadMomentum(fileId)
+   weight = ReadWeight(fileId)
+   time=ReadTime(fileId)
+
+   hyp=sqrt(position[0]**2+position[1]**2+position[2]**2)
+   position /= hyp
+   hyp=sqrt(momentum[0]**2+momentum[1]**2+momentum[2]**2)
+   momentum /= hyp
+   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)
+   else:
+      cond=(abs(costheta)<=1) & (time>0)
+
+   prop = float(shape(time[time<0])[0])/shape(time)[0]
+   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]
+
+def delay_vs_angle_histo(fileId,gen=-1):
+   theta, time, weight = delay_vs_angle(fileId,gen)
+   nbBins = 100
+
+   theta=log10(theta)
+   dt,dtheta=histogram(theta,nbBins,weights=time)
+   dN,dtheta=histogram(theta,nbBins,weights=weight)
+   dtheta=10**dtheta
+   thetacenter=(dtheta[1:nbBins+1]+dtheta[0:nbBins])/2
+   timing = dt/dN
+
+   return thetacenter, timing
+
+def drawDelay_vs_angle(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)
+   
+   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 = delay_vs_angle_histo(fileId)
+   ax.plot(angle,delay/Nmax,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)
+
+   ax.set_xscale('log')
+   ax.set_yscale('log')
+   ax.grid(b=True,which='major')
+   ax.legend(loc="best")
+   ax.set_xlabel("$\\theta $[deg]")
+   ax.set_ylabel("Time delay [s]")
+
+   show()
@@ -4,6 +4,11 @@ from Read import ReadGeneration, ReadWeight, ReadEnergy
  
  
 def drawWeightHisto(fileId):
+   '''
+      Plot particles'weights histogram, generation by generation
+      Input: directory name
+      Output: particles'weights histogram
+   '''
    fig = figure()
    ax = fig.add_subplot(111)
    nbBins =100 
@@ -24,6 +29,11 @@ def drawWeightHisto(fileId):
  
 #######################################################################################
 def weightVsEnergy(fileId,gen=-1):
+   '''
+      Compute particles'weight versus energy
+      Input: directory name
+      Input (optional): generation (default = all)
+   '''
    energy = ReadEnergy(fileId)
    weight = ReadWeight(fileId)
    nbBins = 100 
@@ -43,6 +53,11 @@ def weightVsEnergy(fileId,gen=-1):
    return enercenter, weight
  
 def drawWeightVsEnergy(fileId):
+   '''
+      Plot particles'weight versus energy, generation by generation
+      Input: directory name
+      Output: graph of particles'weight versus energy
+   '''
    fig = figure()
    ax = fig.add_subplot(111)