Correction on the computation of the spectrum for generation==0

Thomas Fitoussi
1 parent 631c85a9
Showing 8 changed files with 63 additions and 31 deletions Show diff stats
Analysis.py
Modules/Angle.pyc
Modules/Generation.py
Modules/Generation.pyc
Modules/Spectrum.py
Modules/Spectrum.pyc
Modules/Timing.py
Modules/Timing.pyc
@@ -7,14 +7,18 @@ 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()
+def InProgress(fileId,i,Nmax):
+   print "Work on ", fileId, ": ", (i*100)/Nmax, "% done"
-if shape(argv)[0] < 2:
-   error("not enough arguments (at least 1)")
+if shape(argv)[0] < 2: 
+   print "not enough arguments (at least 1)"
+   exit()
+#================================================================#
+#  Parameters
+#================================================================#
 folder = "Results/"
+PowerSpectrum=[1,1.5,2,2.5] 
 for fileId in argv[1:]:
    # read files
@@ -23,6 +27,17 @@ for fileId in argv[1:]:
    weightini, generation, theta_arrival, Esource = ReadExtraFile(fileId,[2,3,4,5])
    nbPhotonsEmitted=ReadProfile(fileId,[3]) 
+   lim =100000
+   time=time[:lim]
+   energy=energy[:lim]
+   weightini=weightini[:lim]
+   generation = generation[:lim]
+   theta_arrival = theta_arrival[:lim]
+   Esource=Esource[:lim]
+
+   cond = (generation==0)
+   test=energy[cond]/Esource[cond]
+
    Gen_contrib = []
    Gen_cont = zeros((int(max(generation))+1))
    Source = []
@@ -31,9 +46,15 @@ for fileId in argv[1:]:
    Angle_Energy = []
    Timing = []
-   for powerlaw_index in [1,1.5,2,2.5]: 
+   for powerlaw_index in PowerSpectrum:
       # apply source spectrum
-      weight = weightini*(Esource)**(1-powerlaw_index)
+      weight_source = Esource**(1-powerlaw_index)
+      weight = weightini* weight_source
+
+      print "%e" %(max(weight_source)/min(weight_source))
+
+      #cond = (energy>95) & (energy<105) & (generation==0)
+      #print weightini[cond], weight[cond]
       #=============================================================================#
       # GENERATION CONTRIBUTION
@@ -52,8 +73,9 @@ for fileId in argv[1:]:
       nbBins = 100
       # draw source spectrum
       Es=array(list(set(Esource)))
+      print shape(Es)
       Ws= (Es)**(1-powerlaw_index)
-      Es,Fs = spectrum(Es,Ws,nbBins)
+      Es,Fs = spectrum(Es,Ws,nbBins=nbBins)
       Es=Es[:,newaxis]
       Fs=Fs[:,newaxis]
       if Source==[]:
@@ -62,11 +84,8 @@ for fileId in argv[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)
+      # primary gamma-rays contribution and full spectrum
+      ener,flux,flux_0 = spectrum(energy,weight,generation,nbBins)
       ener=ener[:,newaxis]
       flux=flux[:,newaxis]
       flux_0=flux_0[:,newaxis]
@@ -113,6 +132,8 @@ for fileId in argv[1:]:
       else:
          Timing = append(Timing,dNdt,axis=1)
+      InProgress(fileId,PowerSpectrum.index(powerlaw_index)+1,shape(PowerSpectrum)[0])
+
    savetxt(folder+fileId+"/Generation.txt",Gen_contrib)
    savetxt(folder+fileId+"/Spectrum.txt",Spectrum)
    savetxt(folder+fileId+"/Source_spectrum.txt",Source)
@@ -15,12 +15,11 @@ def drawGeneration(files):
       i=1
       for powerlaw_index in [1,1.5,2,2.5]: 
          generation, ratio = ReadGeneration(fileId,[0,i])
-         ax.plot(generation-0.5,ratio,drawstyle='steps-mid',label="p=%.1f"%powerlaw_index)
+         ax.plot(generation,ratio,drawstyle='steps-mid',label="p=%.1f"%powerlaw_index)
    xtick=arange(9)
    ax.legend(loc='best')
-   ax.set_xticks(xtick+0.5)
-   ax.set_xticklabels(xtick)
+   ax.set_yscale('log')
    ax.set_xlabel("Generation")
    ax.set_ylabel("Ratio of events [%]")
@@ -5,20 +5,32 @@ from Read import ReadSpectrum, ReadSourceSpectrum
 from Constants import degre
 from Analytic import Ethreshold_gg
-def spectrum(energy,weight,nbBins=100):
+def spectrum(energy,weight,generation=[],nbBins=100):
    '''
       Compute flux versus energy
       Input: events energy and weight (optional: number of bins)
       Output: energy, flux (E^2 dN/dE)
    '''
+   Emin=min(energy)
+   Emax=max(energy)
    energy =log10(energy)
-   flux,ener=histogram(energy,nbBins,weights=weight)
+   flux,ener=histogram(energy,nbBins,range=[log10(Emin),log10(Emax)],weights=weight)
    ener = 10**ener
    enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2
    binSize=ener[1:nbBins+1]-ener[0:nbBins]
    flux = enercenter**2 * flux/binSize 
-   flux /= max(flux) #nbPhotonsEmitted
-   return enercenter,flux
+
+   if generation == []:
+      return enercenter,flux
+   else:
+      cond=(generation==0)
+      flux_0,ener_0=histogram(energy[cond],nbBins,range=[log10(Emin),log10(Emax)],weights=weight[cond])
+      ener_0 = 10**ener_0
+      enercenter_0=(ener_0[1:nbBins+1]+ener_0[0:nbBins])/2
+      binSize=ener_0[1:nbBins+1]-ener_0[0:nbBins]
+      flux_0 = enercenter_0**2 * flux_0/binSize 
+      return enercenter,flux,flux_0
+
 def drawSpectrum(files,PlotAnalytic=False):
    '''
@@ -34,19 +46,19 @@ def drawSpectrum(files,PlotAnalytic=False):
    ax2 = fig.add_subplot(gs[1],sharex=ax1)
    for fileId in files:
-      i=1
-      j=1
+      i=1#3
+      j=1#5
       for powerlaw_index in [1,1.5,2,2.5]: 
          # read files
          Es,Fs = ReadSourceSpectrum(fileId,[0,i])
-         p = ax1.plot(Es,Fs,drawstyle='--')
+         p = ax1.plot(Es,Fs,drawstyle='.')
          # draw full spectrum
          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)
+         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='-.')
+         ax1.plot(ener,flux_0,color=p[0].get_color(),linestyle='--',drawstyle='steps-mid')
          contrib=flux_0/flux
-         ax2.plot(ener,contrib,'-',color=p[0].get_color())
+         ax2.plot(ener,contrib,drawstyle='steps-mid',color=p[0].get_color())
          i+=1
          j+=2
@@ -56,9 +56,9 @@ def drawDelay_vs_angle(fileId):
    ax = fig.add_subplot(111)
    nbBins = 100
-   theta      = ReadExtraFile(fileId,[4])
-   energy     = ReadEnergy(fileId)
-   delay      = ReadTime(fileId)
+   theta  = ReadExtraFile(fileId,[4])
+   energy = ReadEnergy(fileId)
+   delay  = ReadTime(fileId)
    for n in arange(0,8,2):
       Emin = 10**(3-n) #GeV
@@ -74,14 +74,14 @@ def drawDelay_vs_angle(fileId):
    dtheta=10**dtheta
    thetacenter=(dtheta[1:nbBins+1]+dtheta[0:nbBins])/2
    dt=dt/dN
-   ax.plot(thetacenter,dt,color='k',linestyle="steps-mid",linewidth=2)
+   ax.plot(thetacenter,dt,color='m',linestyle="steps-mid",linewidth=2)
    thmin = 1.e-8
    thmax = pi/2.
    nth = 1000
    th = thmin*(thmax/thmin)**(arange(nth)/(nth-1.))
    delay_fit = Analytic_delay_vs_theta(th,fileId)
-   ax.plot(th,delay_fit,'--m',linewidth=2)
+   ax.plot(th,delay_fit,'--k',linewidth=2)
    ax.set_xscale('log')
    ax.set_yscale('log')