Analysis.py used to pre-compute spectrum, timing and radial distribution

computing done for power source spectrum (p= 1, 1.5, 2, 2.5) speed-up graph drawing, no need to recompute each time everything

Analysis.py used to pre-compute spectrum, timing and radial distribution
computing done for power source spectrum (p= 1, 1.5, 2, 2.5) speed-up graph drawing, no need to recompute each time everything
Thomas Fitoussi
1 parent 65135318
Showing 11 changed files with 175 additions and 91 deletions Show diff stats
Analysis.py
Modules/Angle.py
Modules/Angle.pyc
Modules/Generation.py
Modules/Generation.pyc
Modules/Read.py
Modules/Read.pyc
Modules/Spectrum.py
Modules/Spectrum.pyc
Modules/Timing.py
Modules/Timing.pyc
 #!/bin/python
 from sys import argv
-from numpy import shape
+from numpy import append,savetxt, shape, array, newaxis, empty
+from Modules.Read import ReadEnergy, ReadTime, ReadExtraFile, ReadProfile
 from Modules.Spectrum import spectrum
+from Modules.Angle import radial, angle_vs_energy
+from Modules.Timing import timing
 def error(error_type):
    print error_type
    exit()
 if shape(argv)[0] < 2:
-   error("not enough arguments (at least 2)")
+   error("not enough arguments (at least 1)")
-# choose what to print
-print "What will be printed?"
-print "  1. Spectrum"
-print "  2. Image"
-print "  3. Radial distribution VS energy"
-print "  4. Timing distribution VS arrival angle"
-selection = raw_input("selection (number): ")
+folder = "Results/"
-if (selection == "1"): # draw spectrum
+for fileId in argv[1:]:
+   # read files
+   time = ReadTime(fileId)
+   energy = ReadEnergy(fileId)
+   weightini, generation, theta_arrival, Esource = ReadExtraFile(fileId,[2,3,4,5])
+   nbPhotonsEmitted=ReadProfile(fileId,[3]) 
+
+   Source = []
+   Spectrum = []
+   Radial = []
+   Angle_Energy = []
+   Timing = []
+
+   for powerlaw_index in [1,1.5,2,2.5]: 
+      # apply source spectrum
+      weight = weightini*(Esource)**(1-powerlaw_index)
+
+      #=============================================================================#
+      #  SPECTRUM (SOURCE AND MEASURED)
+      #=============================================================================#
+      nbBins = 100
+      # draw source spectrum
+      Es=array(list(set(Esource)))
+      Ws= (Es)**(1-powerlaw_index)
+      Es,Fs = spectrum(Es,Ws,nbBins)
+      Es=Es[:,newaxis]
+      Fs=Fs[:,newaxis]
+      if Source==[]:
+         Source = Es
+         Source = append(Source,Fs,axis=1)
+      else: 
+         Source = append(Source,Fs,axis=1)
+
+      # primary gamma-rays contribution
+      cond = (generation==0)
+      ener,flux_0 = spectrum(energy[cond],weight[cond],nbBins)
+      # full spectrum
+      ener,flux = spectrum(energy,weight,nbBins)
+      ener=ener[:,newaxis]
+      flux=flux[:,newaxis]
+      flux_0=flux_0[:,newaxis]
+      if Spectrum==[]:
+         Spectrum = ener
+         Spectrum = append(Spectrum,flux,axis=1)
+         Spectrum = append(Spectrum,flux_0,axis=1)
+      else:
+         Spectrum = append(Spectrum,flux,axis=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
+      delta_t,dNdt = timing(time,weight,nbBins)
+      delta_t=delta_t[:,newaxis]
+      dNdt=dNdt[:,newaxis]
+      if Timing==[]:
+         Timing = delta_t
+         Timing = append(Timing,dNdt,axis=1)
+      else:
+         Timing = append(Timing,dNdt,axis=1)
+
+   savetxt(folder+fileId+"/Spectrum.txt",Spectrum)
+   savetxt(folder+fileId+"/Source_spectrum.txt",Source)
+   savetxt(folder+fileId+"/Angle_vs_Energy.txt",Angle_Energy)
+   savetxt(folder+fileId+"/Radial_distribution.txt",Radial)
+   savetxt(folder+fileId+"/Timing.txt",Timing)
 from numpy import log10, histogram
 from matplotlib.pyplot import figure, show
 from matplotlib import gridspec
-from Read import ReadEnergy, ReadExtraFile, ReadProfile
+from Read import ReadRadial, ReadAngleVsEnergy
 from Analytic import Analytic_theta
 from Constants import degre
@@ -23,6 +23,24 @@ 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
@@ -33,27 +51,11 @@ def drawAngle_vs_energy(files):
    ax1 = fig.add_subplot(111)
    for fileId in files:
-      energy = ReadEnergy(fileId)
-      weightini, theta_arrival, Esource = ReadExtraFile(fileId,[2,4,5])
-      nbPhotonsEmitted=ReadProfile(fileId,[3]) 
-
-      # apply source spectrum
-      powerlaw_index = 1.5
-      weight = weightini*(Esource/min(Esource))**(1-powerlaw_index)
-      ener,angle = angle_vs_energy(theta_arrival,energy,weight,nbPhotonsEmitted)
-      p=ax1.plot(ener,angle,'-',label="p=%.1f"%powerlaw_index)
-
-      # apply source spectrum
-      powerlaw_index = 2
-      weight = weightini*(Esource/min(Esource))**(1-powerlaw_index)
-      ener,angle = angle_vs_energy(theta_arrival,energy,weight,nbPhotonsEmitted)
-      p=ax1.plot(ener,angle,'-',label="p=%.1f"%powerlaw_index)
-
-      # apply source spectrum
-      powerlaw_index = 2.5
-      weight = weightini*(Esource/min(Esource))**(1-powerlaw_index)
-      ener,angle = angle_vs_energy(theta_arrival,energy,weight,nbPhotonsEmitted)
-      p=ax1.plot(ener,angle,'-',label="p=%.1f"%powerlaw_index)
+      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)
       yfit = Analytic_theta(ener,fileId)
       ax1.plot(ener,yfit,'--',color="m",linewidth=2,label="analytic")
@@ -68,42 +70,19 @@ def drawAngle_vs_energy(files):
    show()
-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 drawRadial(files,thetamax=90):
+def drawRadial(files):
    '''
       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:
-      weightini, theta_arrival, Esource = ReadExtraFile(fileId,[2,4,5])
-      nbPhotonsEmitted=ReadProfile(fileId,[3]) 
-
-      for powerlaw_index in [1.5,2,2.5]: 
-         # apply source spectrum
-         weight = weightini*(Esource/min(Esource))**(1-powerlaw_index)
-         theta,dndtheta2=radial(theta_arrival,weight,nbPhotonsEmitted,thetamax)
+      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")
 from numpy import shape, arange
 from matplotlib.pyplot import figure, show, hist
-from Read import ReadGeneration, ReadWeight
+from Read import ReadExtraFile
 def drawGeneration(files):
@@ -13,8 +13,7 @@ def drawGeneration(files):
    ax = fig.add_subplot(111)
    nbBins = 7
    for fileId in files:
-      weight = ReadWeight(fileId)
-      generation = ReadGeneration(fileId)
+      weight, generation = ReadExtraFile(fileId,[2,3])
       print "file", fileId, "->", shape(generation)[0], "events" 
       hist(generation,nbBins,range=[2,9],log=1,alpha=.5, weights=weight,label=fileId)
@@ -26,3 +25,15 @@ def drawGeneration(files):
    ax.set_ylabel("Number of events")
    show()
+
+def GenContribution(fileId):
+   '''
+      compute the contribution to the final events of each
+      generation
+   '''
+   weight, generation = ReadExtraFile(fileId,[2,3])
+   
+   NbTotEvents = sum(weight)
+
+   for gen in list(set(generation)):
+      print "Generation", gen, " -> ", sum(weight[generation==gen])/NbTotEvents *100, "\percent"
@@ -59,3 +59,18 @@ def ReadElmag(fileName):
    energy *= 10**(-9) #GeV
    flux=(fluxGamma+fluxLepton)*1e-9 #GeV normalized to 1 initial photon
    return energy,flux
+
+def ReadSpectrum(fileName,cols=[0,1,2,3,4,5,6,7,8]): 
+   return loadtxt(resultDirectory+fileName+"/Spectrum.txt",unpack=True,usecols=cols)
+
+def ReadSourceSpectrum(fileName,cols=[0,1,2,3,4]): 
+   return loadtxt(resultDirectory+fileName+"/Source_spectrum.txt",unpack=True,usecols=cols)
+
+def ReadAngleVsEnergy(fileName,cols=[0,1,2,3,4]): 
+   return loadtxt(resultDirectory+fileName+"/Angle_vs_Energy.txt",unpack=True,usecols=cols)
+
+def ReadRadial(fileName,cols=[0,1,2,3,4]): 
+   return loadtxt(resultDirectory+fileName+"/Radial_distribution.txt",unpack=True,usecols=cols)
+
+def ReadTiming(fileName,cols=[0,1,2,3,4]): 
+   return loadtxt(resultDirectory+fileName+"/Timing.txt",unpack=True,usecols=cols)
 from numpy import histogram, log10, sqrt, array
 from matplotlib.pyplot import figure, show, setp
 from matplotlib import gridspec
-from Read import ReadEnergy, ReadExtraFile, ReadProfile
+from Read import ReadSpectrum, ReadSourceSpectrum
 from Constants import degre
 from Analytic import Ethreshold_gg
-def spectrum(energy,weight,nbPhotonsEmitted=1,nbBins=100):
+def spectrum(energy,weight,nbBins=100):
    '''
       Compute flux versus energy
       Input: events energy and weight (optional: number of bins)
@@ -17,7 +17,7 @@ def spectrum(energy,weight,nbPhotonsEmitted=1,nbBins=100):
    enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2
    binSize=ener[1:nbBins+1]-ener[0:nbBins]
    flux = enercenter**2 * flux/binSize 
-   flux /= nbPhotonsEmitted
+   flux /= max(flux) #nbPhotonsEmitted
    return enercenter,flux
 def drawSpectrum(files,PlotAnalytic=False):
@@ -34,23 +34,22 @@ def drawSpectrum(files,PlotAnalytic=False):
    ax2 = fig.add_subplot(gs[1],sharex=ax1)
    for fileId in files:
-      # read files
-      energy = ReadEnergy(fileId)
-      weightini, generation, theta_arrival, Esource = ReadExtraFile(fileId,[2,3,4,5])
-      nbPhotonsEmitted=ReadProfile(fileId,[3]) 
-
-      for powerlaw_index in [1.5,2,2.5]: 
-         # apply source spectrum
-         weight = weightini*(Esource/min(Esource))**(1-powerlaw_index)
+      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,drawstyle='--')
          # draw full spectrum
-         ener,flux = spectrum(energy,weight,nbPhotonsEmitted)
-         p = ax1.plot(ener,flux,drawstyle='-',label="p=%.1f"%powerlaw_index)
+         ener,flux,flux_0 = ReadSpectrum(fileId,[0,j,j+1])
+         ax1.plot(ener,flux,color=p[0].get_color(),drawstyle='-',label="p=%.1f"%powerlaw_index)
          # primary gamma-rays contribution
-         cond = (generation==0)
-         ener,flux_cond = spectrum(energy[cond],weight[cond],nbPhotonsEmitted)
-         ax1.plot(ener,flux_cond,color=p[0].get_color(),linestyle='-.')
-         contrib=flux_cond/flux
+         ax1.plot(ener,flux_0,color=p[0].get_color(),linestyle='-.')
+         contrib=flux_0/flux
          ax2.plot(ener,contrib,'-',color=p[0].get_color())
+         
+         i+=1
+         j+=2
       # add spectrum for Fermi/LAT and CTA (FoV?)
 from numpy import histogram, log10, shape, arange, pi
 from matplotlib.pyplot import figure, show
-from Read import ReadTime, ReadExtraFile, ReadProfile, ReadEnergy
+from Read import ReadTime, ReadTiming, ReadExtraFile, ReadProfile, ReadEnergy
 from Constants import yr, degre
 from Analytic import Analytic_delay_vs_theta
-def timing(time,weight,NbPhotonsEmitted=1,nbBins=100):
+def timing(time,weight,nbBins=100):
    '''
       Compute flux versus time delay
       Input: directory name
@@ -23,7 +23,7 @@ def timing(time,weight,NbPhotonsEmitted=1,nbBins=100):
    dN,dt=histogram(time,nbBins,weights=weight)
    timecenter=(dt[1:nbBins+1]+dt[0:nbBins])/2
    binSize=dt[1:nbBins+1]-dt[0:nbBins]
-   dNdt=dN/(NbPhotonsEmitted*binSize)
+   dNdt=dN/(max(dN)*binSize)
    return 10**timecenter, dNdt
@@ -36,14 +36,9 @@ def drawTiming(files):
    fig = figure()
    ax = fig.add_subplot(111)
    for fileId in files:
-      time = ReadTime(fileId)
-      weightini, Esource = ReadExtraFile(fileId,[2,5])
-      nbPhotonsEmitted=ReadProfile(fileId,[3]) 
-
-      for powerlaw_index in [1.5,2,2.5]: 
-         # apply source spectrum
-         weight = weightini*(Esource/min(Esource))**(1-powerlaw_index)
-         delta_t,dNdt = timing(time,weight,nbPhotonsEmitted)
+      i=1
+      for powerlaw_index in [1,1.5,2,2.5]: 
+         delta_t,dNdt = ReadTiming(fileId,[0,i])
          ax.plot(delta_t,dNdt,linestyle="steps-mid",label="p=%.1f"%powerlaw_index)
    ax.set_xscale('log')
@@ -51,7 +46,7 @@ def drawTiming(files):
    ax.grid(b=True,which='major')
    ax.legend(loc="best")
    ax.set_xlabel("Time delay [s]")
-   ax.set_ylabel("$t\ dN/dt$")
+   ax.set_ylabel("$t\ dN/dt$ [arbitrary units]")
    show()