From 0b56ded50cdb37ff0f69daba43d96cd12425405e Mon Sep 17 00:00:00 2001
From: Thomas Fitoussi <tfitoussi@mailoo.org>
Date: Wed, 26 Aug 2015 14:38:56 +0200
Subject: [PATCH] take into account source angle of emission    => selection for a cone emission rather than isotropic

---
 Analysis.py          | 79 +++++++++++++++++++++++++++++++++++++++----------------------------------------
 Modules/Map.py       | 96 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------
 Modules/Map.pyc      | Bin 4925 -> 0 bytes
 Modules/Spectrum.py  | 35 +++++++++++++++++++++++------------
 Modules/Spectrum.pyc | Bin 3395 -> 0 bytes
 5 files changed, 118 insertions(+), 92 deletions(-)

diff --git a/Analysis.py b/Analysis.py
index fb90c7f..3abd2e3 100755
--- a/Analysis.py
+++ b/Analysis.py
@@ -5,7 +5,7 @@ from numpy import append, savetxt, shape, array, newaxis, zeros, arange
 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, isotrop
+from Modules.Map import computeMap
 from Modules.Angle import angle_vs_energy
 from Modules.Timing import timing
 from Modules.Constants import degre
@@ -19,6 +19,7 @@ if shape(argv)[0] < 2:
 
 #=============================================================================#
 PowerSpectrum=[1,1.5,2,2.5] 
+Jet_Opening=[180,30,60] # degre (180 <=> isotrop)
 
 for fileId in argv[1:]:
    print "#=============================================================================#"
@@ -28,14 +29,14 @@ for fileId in argv[1:]:
    print "# 1. Reading data" 
    time = ReadTime(fileId)
    energy = ReadEnergy(fileId)
-   weightini, generation, theta_arrival, Esource = ReadExtraFile(fileId,[2,3,4,5])
+   weightini, generation, theta_arrival, Esource, dir_source = ReadExtraFile(fileId,[2,3,4,5,6])
    nbPhotonsEmitted=ReadProfile(fileId,[3]) 
+   weightini /= nbPhotonsEmitted
 
    Gen_contrib = []
    Gen_cont = zeros((int(max(generation))+1))
    Source = []
    Spectrum = []
-   Angle_Energy = []
    Timing = []
 
    print "# 2. Computing powerlaw spectrum" 
@@ -61,7 +62,7 @@ for fileId in argv[1:]:
       nbBins = 50
       # draw source spectrum
       Es=array(list(set(Esource)))
-      Ws= (Es/min(Es))**(1-powerlaw_index)
+      Ws= (Es/min(Es))**(1-powerlaw_index) /nbPhotonsEmitted
       Es,Fs = spectrum(Es,Ws,nbBins=nbBins)
       Es=Es[:,newaxis]
       Fs=Fs[:,newaxis]
@@ -102,44 +103,42 @@ for fileId in argv[1:]:
    #=============================================================================#
    #  IMAGING, RADIAL DISTRIBUTION AND ANGLE VERSUS ENERGY
    #=============================================================================#
-   print "# 3. compute images and radial distribution" 
-   powerlaw_index = 2
-   Elim = 1e-1 # GeV
-   thetalim = 20 # degre
-   nbBins = 100
+   Angle_Energy = []
    Radial = []
+   print "# 3. compute images and radial distribution" 
+   for jet_opening_angle in Jet_Opening:
+      powerlaw_index = 2
+      Elim = 1e-1 # GeV
+      thetalim = 20 # degre
+      nbBins = 100
+
+      thetaDir,phiDir = ReadMomentumAngle(fileId)*degre
+      thetaPos,phiPos = ReadPositionAngle(fileId)*degre
+      theta = thetaDir - thetaPos
+      phi = phiDir -phiPos
+
+      # apply source spectrum
+      weight_source = (Esource/min(Esource))**(1-powerlaw_index)
+      weight = weightini* weight_source
+
+      # apply selection
+      theta2,dndtheta2,ener,angle  = computeMap(theta,phi,weight,energy,dir_source,fileId,
+            jet_opening_angle,Elim,nbBins,borne=[thetalim,thetalim])
+      theta2=theta2[:,newaxis]
+      dndtheta2=dndtheta2[:,newaxis]
+      ener=ener[:,newaxis]
+      angle=angle[:,newaxis]
+
+      if Radial==[]:
+         Radial = theta2
+         Radial = append(Radial,dndtheta2,axis=1)
+         Angle_Energy = ener
+         Angle_Energy = append(Angle_Energy,angle,axis=1)
+      else:
+         Radial = append(Radial,dndtheta2,axis=1)
+         Angle_Energy = append(Angle_Energy,angle,axis=1)
 
-   thetaDir,phiDir = ReadMomentumAngle(fileId)*degre
-   thetaPos,phiPos = ReadPositionAngle(fileId)*degre
-   theta = thetaDir - thetaPos
-   phi = phiDir -phiPos
-
-   # apply source spectrum
-   weight_source = (Esource/min(Esource))**(1-powerlaw_index)
-   weight = weightini* weight_source
-
-   # apply selection
-   cond=(energy>Elim) #& (abs(theta)<thetalim) & (abs(phi)<thetalim)
-   theta=theta[cond]
-   phi=phi[cond]
-   weight=weight[cond]
-   energy=energy[cond]
-   theta2,dndtheta2  = computeMap(theta,phi,weight,energy,fileId,nbBins,
-         source=isotrop(),borne=[thetalim,thetalim])
-   ener,angle = angle_vs_energy(theta,energy,weight,nbBins)
-   theta2=theta2[:,newaxis]
-   dndtheta2=dndtheta2[:,newaxis]
-   ener=ener[:,newaxis]
-   angle=angle[:,newaxis]
-
-   if Radial==[]:
-      Radial = theta2
-      Radial = append(Radial,dndtheta2,axis=1)
-      Angle_Energy = ener
-      Angle_Energy = append(Angle_Energy,angle,axis=1)
-   else:
-      Radial = append(Radial,dndtheta2,axis=1)
-      Angle_Energy = append(Angle_Energy,angle,axis=1)
+      InProgress(Jet_Opening.index(jet_opening_angle)+1,shape(Jet_Opening)[0])
 
    #=============================================================================#
    print "# 4. writing files" 
diff --git a/Modules/Map.py b/Modules/Map.py
index 5cfff15..8e05485 100644
--- a/Modules/Map.py
+++ b/Modules/Map.py
@@ -1,6 +1,5 @@
 from numpy import shape, histogram2d, arange, newaxis, array, nditer, zeros, histogram
 from numpy import sqrt, pi, exp, log, log10, sin, cos, where
-from scipy.signal import convolve2d
 from matplotlib.pyplot import figure, matshow, savefig, show
 from matplotlib.colors import LogNorm
 from Read import  resultDirectory, ReadRadial, ReadAngleVsEnergy
@@ -15,16 +14,7 @@ philim=float(philim)
 thetalim=float(thetalim)
 limits = [[-thetalim,thetalim],[-philim,philim]]
 
-def isotrop():
-   return array([[1]]), "isotropic"
-
-def Gaussian2D(center=[0,0],fwhm=0.01):
-   x = arange(-thetalim,thetalim+step,step)
-   y = arange(-philim,philim+step,step)
-   y = y[:,newaxis]
-   return exp(-4*log(2) * ((x-center[0])**2 + (y+center[1])**2) / fwhm**2), "jet_"+str(fwhm)
-
-def computeMap(theta,phi,weight,energy,saveId,nbBins=50,source=isotrop(),borne=[180,180]):
+def computeMap(theta,phi,weight,energy,dir_source,saveId,jet_opening_angle=180,Elim=1e-3,nbBins=50,borne=[180,180]):
    '''
       Plot 2D histogram of arrival direction of the photons with energy upper than Elim
       Input: list of directories
@@ -33,20 +23,22 @@ def computeMap(theta,phi,weight,energy,saveId,nbBins=50,source=isotrop(),borne=[
       Input (optional): limits on theta, phi 
       Output: 2D histogram of theta cos(phi), theta sin(phi)
    '''
-   fig = figure()
 
-   H,xedges,yedges = histogram2d(theta*cos(phi),theta*sin(phi),bining,weights=weight,range=limits)
-   conv=convolve2d(H,source[0],boundary="wrap",mode="same")
-   #H2,xedges,yedges = histogram2d(theta*cos(phi),theta*sin(phi),bining,weights=energy,range=limits)
-   #conv2=convolve2d(H2,source[0],boundary="wrap",mode="same")
+   # apply selection
+   cond = (dir_source*degre <= jet_opening_angle) & (energy>Elim) #& (abs(theta)<thetalim) & (abs(phi)<thetalim)
+   theta    = theta[cond]
+   phi      = phi[cond]
+   weight   = weight[cond]
+   energy   = energy[cond]
+
+   fig = figure()
 
    #===============================================================================================#
    #  IMAGE
    #===============================================================================================#
    ax = fig.add_subplot(111)
-   im = ax.matshow(conv,norm=LogNorm(),aspect='auto',
-                     extent=[xedges[0],xedges[-1],yedges[0],yedges[-1]])
-   #H,xedges,yedges,im = ax.hist2d(theta,phi,bining,weights=weight,norm=LogNorm())
+   H,xedges,yedges,im = ax.hist2d(theta*cos(phi),theta*sin(phi),bining,weights=weight,
+         norm=LogNorm(),range=limits)
    cbar=fig.colorbar(im, ax=ax)
    cbar.ax.set_ylabel("counts")
    ax.set_xlim((-borne[0],borne[0]))
@@ -54,30 +46,32 @@ def computeMap(theta,phi,weight,energy,saveId,nbBins=50,source=isotrop(),borne=[
    ax.set_xlabel("$\\theta$ $\\cos(\\phi)$ [deg]")
    ax.set_ylabel("$\\theta$ $\\sin(\\phi)$ [deg]")
    ax.grid(b=True,which='major')
-   ax.set_title(saveId+" - "+source[1])
-
-   savefig(resultDirectory+saveId+"/map_"+source[1]+".png")
+   if jet_opening_angle == 180:
+      ax.set_title(saveId+" - isotrop")
+      savefig(resultDirectory+saveId+"/map_isotrop.png")
+   else:
+      ax.set_title(saveId+" - jet opening angle: %d degre"%(int(jet_opening_angle)))
+      savefig(resultDirectory+saveId+"/map_"+str(int(jet_opening_angle))+".png")
 
+   if borne != [180,180]: 
+      H,xedges,yedges,im = ax.hist2d(theta*cos(phi),theta*sin(phi),bining,weights=weight,
+            norm=LogNorm(),range=limits)
    #===============================================================================================#
    #  RADIAL DISTRIBUTION AND ENERGY VERSUS ANGLE
    #===============================================================================================#
-   theta2 = zeros(shape(conv)[0]*shape(conv)[1])
-   N_evts = zeros(shape(conv)[0]*shape(conv)[1])
-   #E_evts = zeros(shape(conv)[0]*shape(conv)[1])
+   theta2 = zeros(shape(H)[0]*shape(H)[1])
+   N_evts = zeros(shape(H)[0]*shape(H)[1])
    theta_cos_phi = xedges[::-1]
    theta_sin_phi = yedges[::-1]
-   it = nditer(conv,flags=["multi_index"])
-   #it2 = nditer(conv2,flags=["multi_index"])
+   it = nditer(H,flags=["multi_index"])
    k = 0
    while not it.finished:
       i = int(it.multi_index[1])
       j = int(it.multi_index[0])
       theta2[k]= sqrt(theta_cos_phi[i]**2+theta_sin_phi[j]**2)
       N_evts[k]= it[0]
-      #E_evts[k]= it2[0]
       k += 1
       it.iternext()
-      #it2.iternext()
 
    theta2=log10(theta2)
    thetamax=max(borne)
@@ -87,18 +81,14 @@ def computeMap(theta,phi,weight,energy,saveId,nbBins=50,source=isotrop(),borne=[
    dtheta = angle[1:nbBins+1]-angle[0:nbBins]
    dndtheta = dn/dtheta*thetacenter
 
-   #theta2 = 10**theta2
-   #E_evts =  E_evts[E_evts!=0]
-   #N_evts =  N_evts[N_evts!=0]
-   #theta2 =  theta2[E_evts!=0]
-   #E_evts =log10(E_evts)
-   #angle,ener =histogram(E_evts,nbBins,weights=N_evts*theta2)
-   #nb,ener =histogram(E_evts,nbBins,weights=N_evts)
-   #ener = 10**ener
-   #enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2
-   #angle /= nb*degre 
+   energy =log10(energy)
+   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 /= nb 
 
-   return thetacenter, dndtheta #, enercenter, angle
+   return thetacenter, dndtheta , enercenter, angle
 
 def drawRadial(files):
    '''
@@ -121,3 +111,29 @@ def drawRadial(files):
    ax.set_ylabel("$\\theta*dN/d\\theta$")
    
    show()
+
+def drawAngle_vs_energy(files):
+   '''
+      Plot arrival angle versus energy + analytic expression
+      Input: list of directories
+      Output: graph of arrival angle versus energy
+   '''
+   fig = figure()
+   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)
+
+      yfit = Analytic_theta(ener,fileId)
+      ax1.plot(ener,yfit,'--',color="m",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_xlabel("energy [GeV]")
+
+   show()
diff --git a/Modules/Map.pyc b/Modules/Map.pyc
index 751418e..2d6dfa4 100644
Binary files a/Modules/Map.pyc and b/Modules/Map.pyc differ
diff --git a/Modules/Spectrum.py b/Modules/Spectrum.py
index 12e2cf3..10d0b4a 100644
--- a/Modules/Spectrum.py
+++ b/Modules/Spectrum.py
@@ -32,7 +32,7 @@ def spectrum(energy,weight,generation=[],nbBins=100):
       return enercenter,flux,flux_0
 
 
-def drawSpectrum(files,PlotAnalytic=False):
+def drawSpectrum(files,all_source_spectrum=False,PlotAnalytic=False):
    '''
       Plot flux versus energy
       Input: list of directories
@@ -46,25 +46,36 @@ def drawSpectrum(files,PlotAnalytic=False):
    ax2 = fig.add_subplot(gs[1],sharex=ax1)
 
    for fileId in files:
-      i=1#3
-      j=1#5
-      for powerlaw_index in [1,1.5,2,2.5]: 
+      if all_source_spectrum:
+         i=1
+         j=1
+         for powerlaw_index in [1,1.5,2,2.5]: 
+            # read files
+            Es,Fs = ReadSourceSpectrum(fileId,[0,i])
+            p = ax1.plot(Es,Fs,linestyle=':')
+            # draw full spectrum
+            ener,flux,flux_0 = ReadSpectrum(fileId,[0,j,j+1])
+            ax1.plot(ener,flux,color=p[0].get_color(),drawstyle='steps-mid',label="p=%.1f"%powerlaw_index)
+            # primary gamma-rays contribution
+            ax1.plot(ener,flux_0,color=p[0].get_color(),linestyle='--',drawstyle='steps-mid')
+            contrib=flux_0/flux
+            ax2.plot(ener,contrib,drawstyle='steps-mid',color=p[0].get_color())
+            
+            i+=1
+            j+=2
+
+      else:
          # read files
-         Es,Fs = ReadSourceSpectrum(fileId,[0,i])
+         Es,Fs = ReadSourceSpectrum(fileId,[0,3])
          p = ax1.plot(Es,Fs,linestyle=':')
          # draw full spectrum
-         ener,flux,flux_0 = ReadSpectrum(fileId,[0,j,j+1])
-         ax1.plot(ener,flux,color=p[0].get_color(),drawstyle='steps-mid',label="p=%.1f"%powerlaw_index)
+         ener,flux,flux_0 = ReadSpectrum(fileId,[0,5,6])
+         ax1.plot(ener,flux,color=p[0].get_color(),drawstyle='steps-mid',label=fileId)
          # primary gamma-rays contribution
          ax1.plot(ener,flux_0,color=p[0].get_color(),linestyle='--',drawstyle='steps-mid')
          contrib=flux_0/flux
          ax2.plot(ener,contrib,drawstyle='steps-mid',color=p[0].get_color())
          
-         i+=1
-         j+=2
-
-      # add spectrum for Fermi/LAT and CTA (FoV?)
-
    if PlotAnalytic:
       minEner=min(ener)
       alpha=flux[ener==minEner]/sqrt(minEner)
diff --git a/Modules/Spectrum.pyc b/Modules/Spectrum.pyc
index 02e040f..094b483 100644
Binary files a/Modules/Spectrum.pyc and b/Modules/Spectrum.pyc differ
--
libgit2 0.21.2