Commit 374f3d52ed4a3587b3d9f2cb7c21bbd0af4787a1

Authored by Thomas Fitoussi
1 parent 1db105d4
Exists in master

add Weight distribution and versus energy

spectrum can be used with generation selection
Analysis.py
... ... @@ -2,7 +2,7 @@
2 2  
3 3 from sys import argv
4 4 from numpy import shape
5   -from Modules.Spectrum import drawSpectrum
  5 +from Modules.Spectrum import drawSpectrum, drawSpectrum_normaLepton
6 6 from Modules.Angle import drawAngle, drawRadial
7 7 from Modules.Timing import drawTiming
8 8 from Modules.Generation import drawGeneration
... ... @@ -25,9 +25,12 @@ def error(error_type):
25 25 if shape(argv)[0] < 3:
26 26 error("not enough argument")
27 27  
28   -if argv[1] == "spectrum":
  28 +elif argv[1] == "spectrum":
29 29 drawSpectrum(argv[2:])
30 30  
  31 +elif argv[1] == "spectrumLept":
  32 + drawSpectrum_normaLepton(argv[2:])
  33 +
31 34 elif argv[1] == "angle":
32 35 drawAngle(argv[2:])
33 36  
... ...
EBL-spectrums.py
... ... @@ -3,7 +3,7 @@
3 3 import numpy as np
4 4 from scipy.optimize import curve_fit
5 5 import matplotlib.pyplot as plt
6   -from Constantes import *
  6 +from Modules.Constants import *
7 7  
8 8 fig = plt.figure()
9 9 ax = fig.add_subplot(111)
... ...
Modules/Analytic.py
... ... @@ -22,8 +22,8 @@ def Analytic_theta(E_ic,fileId):
22 22 return abs(arcsin(delta_to_theta*sin(delta))*degre)
23 23  
24 24 # Compton accumulation
25   -def ECompton_threshold(Compton_threshold = 0.05):
26   - return Compton_threshold/(4/3*2.7*Tcmbp) *me*1e-6 #GeV
  25 +def ECompton_threshold(Compton_threshold = 0.005):
  26 + return Compton_threshold/(4/3*Ecmb/me*1e-3) *me*1e-6 #GeV
27 27  
28 28 # Compton scattering
29 29 def Dic(Ee=Ee_default): # Ee (GeV)
... ...
Modules/Analytic.pyc
No preview for this file type
Modules/Read.py
... ... @@ -40,8 +40,8 @@ def ReadGeneration(fileName):
40 40 def ReadElmag(fileName):
41 41 Elmag=resultDirectory+fileName+ElmagFile
42 42 energy,fluxGamma,fluxLepton = loadtxt(Elmag,unpack=True,usecols=[0,1,2])
43   - energy *= 10**(-9)
44   - flux=fluxGamma+fluxLepton
  43 + energy *= 10**(-9) #GeV
  44 + flux=(fluxGamma+fluxLepton)*1e-9 #GeV normalized to 1 initial photon
45 45 return energy,flux
46 46  
47 47 def ReadE_source(fileName): # Emax (GeV), redshift, distance (Mpc), Nb photons
... ...
Modules/Read.pyc
No preview for this file type
Modules/Spectrum.py
1   -from numpy import histogram, log10, shape, logspace
  1 +from numpy import histogram, log10, shape, logspace, sqrt
2 2 from matplotlib.pyplot import figure, show, ylim, setp
3 3 from matplotlib import gridspec
4 4 from os.path import isfile
5 5 from Read import ReadEnergy, ReadWeight, ReadElmag, ReadNbLeptonsProd, resultDirectory, ElmagFile
6   -from Analytic import Analytic_flux
  6 +from Read import ReadNphot_source, ReadGeneration
  7 +from Analytic import Analytic_flux, Ethreshold_gg
  8 +
  9 +def spectrum(fileId,gen=-1):
  10 + energy = ReadEnergy(fileId)
  11 + weight = ReadWeight(fileId)
  12 +
  13 + if gen != -1:
  14 + generation = ReadGeneration(fileId)
  15 + energy=energy[generation==gen]
  16 + weight=weight[generation==gen]
  17 + print "file", fileId, "gen=", gen,"->", shape(energy)[0], "events"
  18 + else:
  19 + print "file", fileId, "->", shape(energy)[0], "events"
7 20  
8   -def spectrum(fileName):
9   - energy = ReadEnergy(fileName)
10   - weight = ReadWeight(fileName)
11   - print "file", fileName, "->", shape(energy)[0], "events"
12 21 nbBins=100
13 22 energy =log10(energy)
14 23 flux,ener=histogram(energy,nbBins,weights=weight)
... ... @@ -16,9 +25,46 @@ def spectrum(fileName):
16 25 enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2
17 26 binSize=ener[1:nbBins+1]-ener[0:nbBins]
18 27 flux = enercenter**2 * flux/binSize
  28 + maxflux = ReadNphot_source(fileId)
  29 + flux /= maxflux
19 30 return enercenter,flux
20 31  
21   -def drawSpectrum(files):
  32 +###########################################################################################
  33 +def drawSpectrum(files,PlotElmag=False,PlotAnalytic=False):
  34 + fig = figure()
  35 + ax = fig.add_subplot(111)
  36 +
  37 + for fileId in files:
  38 + energy,flux = spectrum(fileId)
  39 + if PlotElmag:
  40 + maxflux = max(flux)*(.8)
  41 + p=ax.plot(energy,flux,drawstyle='steps-mid',label=fileId)
  42 +
  43 + if PlotElmag & isfile(resultDirectory+fileId+ElmagFile):
  44 + energy, flux = ReadElmag(fileId)
  45 + maxflux = max(flux)
  46 + flux /= maxflux
  47 + ax.plot(energy,flux,"--",color=p[0].get_color(),label="Elmag - "+fileId)
  48 +
  49 + if PlotAnalytic:
  50 + alpha=5e3
  51 + ax.plot(energy,alpha*sqrt(energy),'--m',linewidth=2)
  52 + ax.plot(energy,energy**(0)*max(flux),'--m',linewidth=2)
  53 + ax.axvline(x=Ethreshold_gg(), ymin=0., ymax = 1., color='r', linewidth=2)
  54 +
  55 + ax.set_xscale('log')
  56 + ax.set_yscale('log')
  57 + ax.grid(b=True,which='major')
  58 + ax.legend(loc="best")
  59 + if PlotElmag:
  60 + ax.set_ylabel("$E^2 \\frac{dN}{dE}$ [normalized]")
  61 + else:
  62 + ax.set_ylabel("$E^2 \\frac{dN}{dE}$ [GeV / photon]")
  63 + ax.set_xlabel("energy [GeV]")
  64 +
  65 + show()
  66 +
  67 +def drawSpectrum_normaLepton(files):
22 68 fig = figure()
23 69 gs = gridspec.GridSpec(2, 1, height_ratios=[4,1])
24 70 fig.subplots_adjust(hspace=0.001)
... ... @@ -37,16 +83,11 @@ def drawSpectrum(files):
37 83 error=flux/yfit-1
38 84 ax2.plot(energy,error,'+',color=p[0].get_color())
39 85  
40   - Elmag=resultDirectory+fileId+ElmagFile
41   - if isfile(Elmag):
42   - energy, flux = ReadElmag(fileId)
43   - ax1.plot(energy,flux,"-",color=p[0].get_color(),label="Elmag - "+fileId)
44   -
45 86 ax1.set_xscale('log')
46 87 ax1.set_yscale('log')
47 88 ax1.grid(b=True,which='major')
48 89 ax1.legend(loc="best")
49   - ax1.set_ylabel("$E^2 \\frac{dN}{dE}$ [GeV]")
  90 + ax1.set_ylabel("$E^2 \\frac{dN}{dE}$ [GeV] / lepton")
50 91 ax2.set_xscale('log')
51 92 ax2.grid(b=True,which='major')
52 93 ax2.set_xlabel("energy [GeV]")
... ... @@ -57,3 +98,30 @@ def drawSpectrum(files):
57 98 setp(xticklabels, visible=False)
58 99  
59 100 show()
  101 +
  102 +########################################################################################
  103 +def drawSpectrumGen(fileId):
  104 + fig = figure()
  105 + ax = fig.add_subplot(111)
  106 +
  107 + for gen in list(set(ReadGeneration(fileId))):
  108 + energy,flux = spectrum(fileId,gen)
  109 + p=ax.plot(energy,flux,drawstyle='steps-mid',label="gen="+str(int(gen)))
  110 +
  111 + energy,flux = spectrum(fileId)
  112 + p=ax.plot(energy,flux,drawstyle='steps-mid',color='k',label="gen = all")
  113 + alpha=2e3
  114 +
  115 + ax.plot(energy,sqrt(energy)*alpha,'--m',linewidth=2)
  116 + ax.plot(energy,energy**(0)*max(flux),'--m',linewidth=2)
  117 + ax.axvline(x=Ethreshold_gg(), ymin=0., ymax = 1., color='r', linewidth=2)
  118 +
  119 + ax.set_xscale('log')
  120 + ax.set_yscale('log')
  121 + ax.grid(b=True,which='major')
  122 + ax.legend(loc="best")
  123 + ax.set_title(fileId)
  124 + ax.set_ylabel("$E^2 \\frac{dN}{dE}$ [GeV / photon]")
  125 + ax.set_xlabel("energy [GeV]")
  126 +
  127 + show()
... ...
Modules/Spectrum.pyc
No preview for this file type
Modules/Weight.py 0 โ†’ 100644
... ... @@ -0,0 +1,66 @@
  1 +from numpy import shape, arange, log10, histogram
  2 +from matplotlib.pyplot import figure, show, hist
  3 +from Read import ReadGeneration, ReadWeight, ReadEnergy
  4 +
  5 +
  6 +def drawWeightHisto(fileId):
  7 + fig = figure()
  8 + ax = fig.add_subplot(111)
  9 + nbBins =100
  10 +
  11 + weight = ReadWeight(fileId)
  12 + generation = ReadGeneration(fileId)
  13 + print "file", fileId, "->", shape(generation)[0], "events"
  14 +
  15 + for gen in list(set(generation)):
  16 + hist(log10(weight[generation==gen]),nbBins,log=1,range=[0,7],alpha=.5,
  17 + label="gen="+str(int(gen)))
  18 +
  19 + ax.legend(loc='best')
  20 + ax.set_xlabel("Weight [log]")
  21 + ax.set_ylabel("Number of events")
  22 +
  23 + show()
  24 +
  25 +#######################################################################################
  26 +def weightVsEnergy(fileId,gen=-1):
  27 + energy = ReadEnergy(fileId)
  28 + weight = ReadWeight(fileId)
  29 + nbBins = 100
  30 +
  31 + if gen != -1:
  32 + generation = ReadGeneration(fileId)
  33 + energy=energy[generation==gen]
  34 + weight=weight[generation==gen]
  35 + nbBins =40
  36 +
  37 + weight,ener=histogram(log10(energy),nbBins,weights=log10(weight))
  38 + nb,ener=histogram(log10(energy),nbBins)
  39 + ener=10**ener
  40 + enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2
  41 + weight /= nb
  42 +
  43 + return enercenter, weight
  44 +
  45 +def drawWeightVsEnergy(fileId):
  46 + fig = figure()
  47 + ax = fig.add_subplot(111)
  48 +
  49 + generation = ReadGeneration(fileId)
  50 + print "file", fileId, "->", shape(generation)[0], "events"
  51 +
  52 + for gen in list(set(generation)):
  53 + energy,weight = weightVsEnergy(fileId,gen)
  54 + ax.plot(energy,weight,".",label="gen="+str(int(gen)))
  55 +
  56 + energy,weight = weightVsEnergy(fileId)
  57 + ax.plot(energy,weight,"k",label="gen=all")
  58 +
  59 + ax.set_xscale('log')
  60 + ax.grid(b=True,which='major')
  61 + ax.legend(loc="best")
  62 + ax.set_title(fileId)
  63 + ax.set_xlabel("energy [GeV]")
  64 + ax.set_ylabel("mean weight [log]")
  65 +
  66 + show()
... ...