Put radial distribution inside the map computation

=> can tke into account the jet opening if so

Put radial distribution inside the map computation
=> can tke into account the jet opening if so
Thomas Fitoussi
1 parent 49a0fe94
Showing 7 changed files with 163 additions and 118 deletions Show diff stats
Analysis.py
Modules/Angle.py
Modules/Angle.pyc
Modules/Generation.py
Modules/Map.py
Modules/Map.pyc
README.md
@@ -5,60 +5,40 @@ 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 printMap, isotrop
-from Modules.Angle import radial, angle_vs_energy
+from Modules.Map import computeMap, isotrop
+from Modules.Angle import angle_vs_energy
 from Modules.Timing import timing
 from Modules.Constants import degre
-def InProgress(fileId,step,i,Nmax):
-   print "Work on ", fileId, "(step "+ step +": ", (i*100)/Nmax, "% done"
+def InProgress(i,Nmax):
+   print "  ", (i*100)/Nmax, "% done"
 if shape(argv)[0] < 2: 
    print "not enough arguments (at least 1)"
    exit()
-#================================================================#
-#  Parameters
-#================================================================#
+#=============================================================================#
 PowerSpectrum=[1,1.5,2,2.5] 
-powerlaw_index = 2
-Elim = 1e-1 # GeV
-thetalim = 20 # degre
 for fileId in argv[1:]:
+   print "#=============================================================================#"
+   print "#    Analysis of", fileId
+   print "#=============================================================================#"
    # read files
+   print "# 1. Reading data" 
    time = ReadTime(fileId)
    energy = ReadEnergy(fileId)
    weightini, generation, theta_arrival, Esource = ReadExtraFile(fileId,[2,3,4,5])
    nbPhotonsEmitted=ReadProfile(fileId,[3]) 
-   #=============================================================================#
-   #  IMAGING 
-   #=============================================================================#
-   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)
-   weight=weight[cond]
-   theta=theta[cond]
-   phi=phi[cond]
-   printMap(theta,phi,weight,fileId,source=isotrop(),borne=[thetalim,thetalim])
-
    Gen_contrib = []
    Gen_cont = zeros((int(max(generation))+1))
    Source = []
    Spectrum = []
-   Radial = []
    Angle_Energy = []
    Timing = []
+   print "# 2. Computing powerlaw spectrum" 
    for powerlaw_index in PowerSpectrum:
       # apply source spectrum
       weight_source = (Esource/min(Esource))**(1-powerlaw_index)
@@ -78,7 +58,7 @@ for fileId in argv[1:]:
       #=============================================================================#
       #  SPECTRUM (SOURCE AND MEASURED)
       #=============================================================================#
-      nbBins = 100
+      nbBins = 50
       # draw source spectrum
       Es=array(list(set(Esource)))
       Ws= (Es/min(Es))**(1-powerlaw_index)
@@ -105,28 +85,6 @@ for fileId in argv[1:]:
          Spectrum = append(Spectrum,flux_0,axis=1)
       #=============================================================================#
-      #  RADIAL DISTRIBUTION AND ANGLE VERSUS ENERGY
-      #=============================================================================#
-      nbBins = 100
-      ener,angle = angle_vs_energy(theta_arrival,energy,weight,nbPhotonsEmitted)
-      ener=ener[:,newaxis]
-      angle=angle[:,newaxis]
-      if Angle_Energy==[]:
-         Angle_Energy = ener
-         Angle_Energy = append(Angle_Energy,angle,axis=1)
-      else:
-         Angle_Energy = append(Angle_Energy,angle,axis=1)
-
-      theta,dndtheta2 = radial(theta_arrival,weight,nbPhotonsEmitted)
-      theta=theta[:,newaxis]
-      dndtheta2=dndtheta2[:,newaxis]
-      if Radial==[]:
-         Radial = theta
-         Radial = append(Radial,dndtheta2,axis=1)
-      else:
-         Radial = append(Radial,dndtheta2,axis=1)
-
-      #=============================================================================#
       #  TIME DISTRIBUTION AND TIME DELAY VERSUS ANGLE
       #=============================================================================#
       nbBins = 100
@@ -139,11 +97,57 @@ for fileId in argv[1:]:
       else:
          Timing = append(Timing,dNdt,axis=1)
-      InProgress(fileId,"1",PowerSpectrum.index(powerlaw_index)+1,shape(PowerSpectrum)[0])
+      InProgress(PowerSpectrum.index(powerlaw_index)+1,shape(PowerSpectrum)[0])
+
+   #=============================================================================#
+   #  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
+   Radial = []
+
+   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)
+
+   #=============================================================================#
+   print "# 4. writing files" 
    savetxt(resultDirectory+fileId+"/Generation.txt",Gen_contrib)
    savetxt(resultDirectory+fileId+"/Spectrum.txt",Spectrum)
    savetxt(resultDirectory+fileId+"/Source_spectrum.txt",Source)
    savetxt(resultDirectory+fileId+"/Angle_vs_Energy.txt",Angle_Energy)
    savetxt(resultDirectory+fileId+"/Radial_distribution.txt",Radial)
    savetxt(resultDirectory+fileId+"/Timing.txt",Timing)
+   print "#=============================================================================#"
+
@@ -23,24 +23,6 @@ def angle_vs_energy(theta,energy,weight,nbPhotonsEmitted=1,nbBins=40):
    return enercenter, angle
-def radial(theta,weight,nbPhotonsEmitted=1,nbBins=50,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)
-   '''
-
-   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*thetacenter
-   dndtheta /= nbPhotonsEmitted
-
-   return thetacenter, dndtheta
-
 def drawAngle_vs_energy(files):
    '''
       Plot arrival angle versus energy + analytic expression
@@ -51,11 +33,8 @@ def drawAngle_vs_energy(files):
    ax1 = fig.add_subplot(111)
    for fileId in files:
-      i=1
-      for powerlaw_index in [1,1.5,2,2.5]: 
-         # apply source spectrum
-         ener,angle = ReadAngleVsEnergy(fileId,[0,i])
-         p=ax1.plot(ener,angle,'-',label="p=%.1f"%powerlaw_index)
+      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")
@@ -69,27 +48,3 @@ def drawAngle_vs_energy(files):
    ax1.set_xlabel("energy [GeV]")
    show()
-
-def drawRadial(files):
-   '''
-      Plot flux versus arrival angle
-      Input: list of directories
-      Output: graph of flux versus angle
-   '''
-   fig = figure()
-   ax = fig.add_subplot(111)
-
-   for fileId in files:
-      i=1
-      for powerlaw_index in [1,1.5,2,2.5]: 
-         theta,dndtheta2 = ReadRadial(fileId,[0,i])
-         ax.plot(theta,dndtheta2,drawstyle='steps-mid',label="p=%.1f"%powerlaw_index)
-
-   ax.legend(loc="best")
-   ax.set_xscale('log')
-   ax.set_yscale('log')
-   ax.grid(b=True,which='major')
-   ax.set_xlabel("$\\theta$ [rad]")
-   ax.set_ylabel("$\\theta*dN/d\\theta$")
-   
-   show()
@@ -2,7 +2,6 @@ from numpy import shape, arange
 from matplotlib.pyplot import figure, show, hist
 from Read import ReadGeneration
-
 def drawGeneration(files):
    '''
       Plot histogram of particles'generations
-from numpy import shape, pi, histogram2d, arange, newaxis, exp, log, array, sin, cos
+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, show, matshow, savefig
+from matplotlib.pyplot import figure, matshow, savefig, show
 from matplotlib.colors import LogNorm
-from Read import  resultDirectory
+from Read import  resultDirectory, ReadRadial, ReadAngleVsEnergy
+from Analytic import Analytic_theta
 from Constants import degre
 thetalim = 180
 philim = 180
 step = 0.25 #deg
-nbBins = array([int(2*philim/step)+1,int(2*thetalim/step)+1])
+bining = array([int(2*philim/step)+1,int(2*thetalim/step)+1])
 philim=float(philim)
 thetalim=float(thetalim)
 limits = [[-thetalim,thetalim],[-philim,philim]]
@@ -22,7 +24,7 @@ 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), "jet_"+str(fwhm)
-def printMap(theta,phi,weight,fileId,source=isotrop(),borne=[180,180]):
+def computeMap(theta,phi,weight,energy,saveId,nbBins=50,source=isotrop(),borne=[180,180]):
    '''
       Plot 2D histogram of arrival direction of the photons with energy upper than Elim
       Input: list of directories
@@ -33,13 +35,18 @@ def printMap(theta,phi,weight,fileId,source=isotrop(),borne=[180,180]):
    '''
    fig = figure()
-   H,xedges,yedges = histogram2d(theta*cos(phi),theta*sin(phi),nbBins,weights=weight,range=limits)
+   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")
+   #===============================================================================================#
+   #  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,nbBins,weights=weight,norm=LogNorm())
+   #H,xedges,yedges,im = ax.hist2d(theta,phi,bining,weights=weight,norm=LogNorm())
    cbar=fig.colorbar(im, ax=ax)
    cbar.ax.set_ylabel("counts")
    ax.set_xlim((-borne[0],borne[0]))
@@ -47,7 +54,70 @@ def printMap(theta,phi,weight,fileId,source=isotrop(),borne=[180,180]):
    ax.set_xlabel("$\\theta$ $\\cos(\\phi)$ [deg]")
    ax.set_ylabel("$\\theta$ $\\sin(\\phi)$ [deg]")
    ax.grid(b=True,which='major')
-   ax.set_title(fileId+" - "+source[1])
+   ax.set_title(saveId+" - "+source[1])
-   savefig(resultDirectory+fileId+"/map_"+source[1]+".png")
-   #show()
+   savefig(resultDirectory+saveId+"/map_"+source[1]+".png")
+
+   #===============================================================================================#
+   #  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])
+   theta_cos_phi = xedges[::-1]
+   theta_sin_phi = yedges[::-1]
+   it = nditer(conv,flags=["multi_index"])
+   #it2 = nditer(conv2,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)
+   dn,angle = histogram(theta2,nbBins,range=[log10(1e-2),log10(thetamax)],weights=N_evts)
+   angle=10**angle
+   thetacenter = (angle[1:nbBins+1]+angle[0:nbBins])/2
+   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 
+
+   return thetacenter, dndtheta #, enercenter, angle
+
+def drawRadial(files):
+   '''
+      Plot flux versus arrival angle
+      Input: list of directories
+      Output: graph of flux versus angle
+   '''
+   fig = figure()
+   ax = fig.add_subplot(111)
+
+   for fileId in files:
+      theta,dndtheta2 = ReadRadial(fileId,[0,1])
+      ax.plot(theta,dndtheta2,drawstyle='steps-mid',label=fileId)
+
+   ax.legend(loc="best")
+   ax.set_xscale('log')
+   ax.set_yscale('log')
+   ax.grid(b=True,which='major')
+   ax.set_xlabel("$\\theta$ [deg]")
+   ax.set_ylabel("$\\theta*dN/d\\theta$")
+   
+   show()
@@ -9,16 +9,31 @@ They also requires:
 Instruction to install it could be found here: http://www.scipy.org/install.html
-## Principle
+## Important remarks
-There is interdependancy between modules. Don't copy one and not the others otherwise they will be
+There is interdependency between modules. Don't copy one and not the others otherwise they will be
 errors.
-In the same time, the function in each module are written in way being used independantly, just
+In the same time, the function in each module are written in way being used independently, just
 import it and python will do the rest. For example, to use the spectrum, use just add to your
 script: "from Modules.Spectrum import spectrum".
-Few modules are available:
+## how-to
+
+* First create a directory call Results in which will be copy the output directory of the
+  simulation. You can named the output directory as you wish. It aims to have separately different
+  simulations and plot them together.
+
+  mkdir Results
+
+* Since data files can become really huge, the "Analysis.py" script aims to pre-generate the need
+  files which will be used later to draw figures.
+
+  python Analysis output
+  or if you want to do more than directory at a time
+  python Analysis output another_result_directory etc
+
+* From you draw some figures using the modules available 
 | Modules         | functions                                                 |
 | --------------- | --------------------------------------------------------- |
@@ -29,8 +44,10 @@ Few modules are available:
 | "Spectrum.py"   | extract and plot the spectrum                             |
 | "timing.py"     | extract and plot delta_t/t                                |
 | "Angle.py"      | extract and plot theta^2 in histogram or versus energy    |
-| "Map.py"        | plot map of the arrival photon, points or histo2d         |
-| "Generation.py" | plot historgram of the arrival particles generation       |
+| "Map.py"        | plot map of the arrival photons histo2d                   |
+| "Generation.py" | plot histogram of the arrival particles generation        |
+
+or develop yours!
 **N.B.:** "draw..." function call a table of files because it allows to plots more than one figure
 on the same graph, rather than the rest of the function can call only one file at a time.