Blame view

Modules/Spectrum.py 2.84 KB
65135318   Thomas Fitoussi   Reset modules
1
2
from numpy import histogram, log10, sqrt, array
from matplotlib.pyplot import figure, show, setp
8ffbee93   Thomas Fitoussi   Analytic expressi...
3
from matplotlib import gridspec
f5a75019   Thomas Fitoussi   Rewrite Analysis....
4
from Read import ReadSpectrum, ReadSourceSpectrum, ReadGeneration
65135318   Thomas Fitoussi   Reset modules
5
from Constants import degre
f5a75019   Thomas Fitoussi   Rewrite Analysis....
6
from Analytic import Ethreshold_gg, Analytic_flux, Eic
374f3d52   Thomas Fitoussi   add Weight distri...
7

b1415bc6   Thomas Fitoussi   Correction of spe...
8
def spectrum(energy,weight,nbBins=100,E_range=[1e-3,1e5]):
f4246390   Thomas Fitoussi   Improve angular d...
9
10
   '''
      Compute flux versus energy
65135318   Thomas Fitoussi   Reset modules
11
12
      Input: events energy and weight (optional: number of bins)
      Output: energy, flux (E^2 dN/dE)
f4246390   Thomas Fitoussi   Improve angular d...
13
   '''
b1415bc6   Thomas Fitoussi   Correction of spe...
14
15
   Emin=E_range[0]#min(energy)
   Emax=E_range[1]#max(energy)
7030f150   Thomas Fitoussi   Full reorganisati...
16
   energy =log10(energy)
b1415bc6   Thomas Fitoussi   Correction of spe...
17
   flux,ener=histogram(energy,nbBins,range=log10(E_range),weights=weight)
7030f150   Thomas Fitoussi   Full reorganisati...
18
19
20
21
   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 
f5a75019   Thomas Fitoussi   Rewrite Analysis....
22
   return enercenter,flux
6633967c   Thomas Fitoussi   Correction on the...
23

f5a75019   Thomas Fitoussi   Rewrite Analysis....
24
def drawSpectrum(files,plot_gen=False,plot_source=False,plot_analytic=False):
f4246390   Thomas Fitoussi   Improve angular d...
25
26
27
   '''
      Plot flux versus energy
      Input: list of directories
65135318   Thomas Fitoussi   Reset modules
28
      Input (optional, default=False): plot analytic expression
f4246390   Thomas Fitoussi   Improve angular d...
29
30
      Output: graph spectrum normalize to 1 source photon
   '''
374f3d52   Thomas Fitoussi   add Weight distri...
31
   fig = figure()
a44230b7   Thomas Fitoussi   Improve approxima...
32
33
34
35
36
   #gs = gridspec.GridSpec(2, 1, height_ratios=[3,1]) 
   #fig.subplots_adjust(hspace=0.001)
   #ax1 = fig.add_subplot(gs[0])
   #ax2 = fig.add_subplot(gs[1],sharex=ax1)
   ax1 = fig.add_subplot(111)
8ffbee93   Thomas Fitoussi   Analytic expressi...
37

7030f150   Thomas Fitoussi   Full reorganisati...
38
   for fileId in files:
f5a75019   Thomas Fitoussi   Rewrite Analysis....
39
40
41
42
43
44
45
46
      # draw full spectrum
      ener,flux = ReadSpectrum(fileId,[0,1])
      p = ax1.plot(ener,flux,drawstyle='steps-mid')

      if plot_gen:
         generation = ReadGeneration(fileId,[0])
         j=2
         for gen in generation:
138898c8   Thomas Fitoussi   Computation of ra...
47
            # read files
f5a75019   Thomas Fitoussi   Rewrite Analysis....
48
49
50
            ener,flux = ReadSpectrum(fileId,[0,j])
            ax1.plot(ener,flux,drawstyle='steps-mid',label="gen=%.0f"%gen)
            j+=1
0b56ded5   Thomas Fitoussi   take into account...
51

f5a75019   Thomas Fitoussi   Rewrite Analysis....
52
53
54
      if plot_source:
         Es,Fs = ReadSourceSpectrum(fileId,[0,1])
         ax1.plot(Es,Fs,color=p[0].get_color(),linestyle=':')
074c1729   Thomas Fitoussi   Analysis.py used ...
55
         
f5a75019   Thomas Fitoussi   Rewrite Analysis....
56
   if plot_analytic:
65135318   Thomas Fitoussi   Reset modules
57
      minEner=min(ener)
f5a75019   Thomas Fitoussi   Rewrite Analysis....
58
59
60
61
62
63
64
65
66
67
68
      #alpha=flux[ener==minEner]/sqrt(minEner)
      #ax1.plot(ener,alpha*sqrt(ener/100),'.-g',linewidth=2)
      ax1.plot(ener,2*Analytic_flux(ener),'--g',linewidth=1)
      Elim= Eic(Ethreshold_gg()/2)
      alpha = 2*4.63e3/Elim**(1/4.)
      ax1.plot(ener,alpha*sqrt(ener),'--r',linewidth=1)
      alpha = 2*14.6e3
      ax1.plot(ener,alpha*(ener/100)**(1/4.),'--m',linewidth=1)
      #ax1.plot(ener,ener**(0)*max(flux)*0.95,'--m',linewidth=2)
      #ax1.axvline(x=Ethreshold_gg(), ymin=0., ymax = 1., color='r', linewidth=2)
      #ax1.axvline(x=1e4, ymin=0., ymax = 1., color='g', linewidth=2)
7030f150   Thomas Fitoussi   Full reorganisati...
69

65135318   Thomas Fitoussi   Reset modules
70
   ax1.legend(loc="best")
8ffbee93   Thomas Fitoussi   Analytic expressi...
71
72
73
   ax1.set_xscale('log')
   ax1.set_yscale('log')
   ax1.grid(b=True,which='major')
65135318   Thomas Fitoussi   Reset modules
74
   ax1.set_ylabel("$E^2 dN/dE$ [GeV]")
a44230b7   Thomas Fitoussi   Improve approxima...
75
76
77
78
79
80
   ax1.set_xlabel("energy [GeV]")
   #ax2.set_xscale('log')
   #ax2.grid(b=True,which='major')
   #ax2.set_xlabel("energy [GeV]")
   #xticklabels = ax1.get_xticklabels()
   #setp(xticklabels, visible=False)
7030f150   Thomas Fitoussi   Full reorganisati...
81
82

   show()