Thomas Fitoussi / simulation-analysis

Commit 7030f150ffda3dcf08de5459d2f03003e170ced5

Authored by Thomas Fitoussi
1 parent e4ba931e
Exists in master

Full reorganisation of the scripts 

writen in modules
   each function can be use independantly rather than full complicate scripts
   /!\ interdependancy between modules
Showing 27 changed files with 412 additions and 447 deletions   Show diff stats
Inline Side-by-side
Analysis.py 0 → 100755
  Diff comments   View file @7030f15
... ... @@ -0,0 +1,37 @@
  1 +#!/bin/python
  2 +
  3 +from sys import argv
  4 +from numpy import shape
  5 +from Modules.Spectrum import drawSpectrum
  6 +from Modules.Angle import drawAngle, drawRadial
  7 +from Modules.Timing import drawTiming
  8 +from Modules.Generation import drawGeneration
  9 +from Modules.Map import drawMap, drawMapHisto
  10 +
  11 +if shape(argv)[0] < 3:
  12 + print "Give at least 2 arguments"
  13 + exit()
  14 +
  15 +if argv[1] == "spectrum":
  16 + drawSpectrum(argv[2:])
  17 +
  18 +elif argv[1] == "angle":
  19 + drawAngle(argv[2:])
  20 +
  21 +elif argv[1] == "radial":
  22 + drawRadial(argv[2:])
  23 +
  24 +elif argv[1] == "timing":
  25 + drawTiming(argv[2:])
  26 +
  27 +elif argv[1] == "generation":
  28 + drawGeneration(argv[2:])
  29 +
  30 +elif argv[1] == "maphisto":
  31 + drawMapHisto(argv[2:])
  32 +
  33 +elif argv[1] == "map":
  34 + drawMap(argv[2:])
  35 +
  36 +else:
  37 + print "bad 1st argument"
... ...
Angle.py deleted
View file @e4ba931
... ... @@ -1,135 +0,0 @@
1   -#!/usr/bin/python
2   -
3   -from sys import argv
4   -import matplotlib.pyplot as plt
5   -from matplotlib import gridspec
6   -from Constantes import *
7   -
8   -if argv[1] == "EGMF":
9   - fileName="Results_EGMF"
10   -elif argv[1] == "EBL":
11   - fileName="Results_EBL"
12   - B=10**(-15)
13   - Model=["Kneiske and Doll - 'best fit'","Kneiske and Doll - 'lower limit'",
14   - "Fraceschini","Finke and Al","Gilmore and Al","Dominguez and Al"]
15   -else:
16   - print "Used: ./Angle_distribution.py arg1 ind1 ind2 ..."
17   - print " arg1 = EGMF or EBL "
18   - print " ind1 = file number \n"
19   - print "Examples: "
20   - print " to study EGMF14: ./Angle_distribution.py EGMF 14 "
21   - print " to study EBL1 and EBL2: ./Angle_distribution.py EBL 1 2 \n"
22   - exit()
23   -
24   -nbBins = 100
25   -Eliminf = 1e0 #GeV
26   -Elimsup = 1e5 #GeV
27   -
28   -Esource=100 #TeV
29   -Dsource=135 #Mpc
30   -delta_to_theta=(lambda_gg/Esource*1e3)/Dsource
31   -
32   -def fit_theta(E,B):
33   - RL0=RL*(Esource/2)*(1e-17/B)
34   - Dic0=Dic/(Esource/2)
35   - delta=Dic0/(2*RL0)*((Esource/2)**2 *Eic/E-1)
36   - return abs(arcsin(delta_to_theta*sin(delta))*degre)
37   -
38   -fig2 = plt.figure()
39   -ax2 = fig2.add_subplot(111)
40   -
41   -fig1 = plt.figure()
42   -gs = gridspec.GridSpec(2, 1, height_ratios=[4,1])
43   -fig1.subplots_adjust(hspace=0.001)
44   -ax11 = fig1.add_subplot(gs[0])
45   -ax12 = fig1.add_subplot(gs[1],sharex=ax11)
46   -
47   -color=['b','r','g','c','m','y']
48   -
49   -ind = 0
50   -for fileId in argv[2:]:
51   - energy,dirx,diry,dirz = loadtxt(fileName+fileId+"/Results_momentum",unpack=True,usecols=[0,1,2,3])
52   - posx,posy,posz = loadtxt(fileName+fileId+"/Results_position",unpack=True,usecols=[1,2,3])
53   - charge,weight,gen=loadtxt(fileName+fileId+"/Results_extra",unpack=True,usecols=[0,2,3])
54   -
55   - hyp=sqrt(posx**2+posy**2+posz**2)
56   - posx=posx/hyp
57   - posy=posy/hyp
58   - posz=posz/hyp
59   - hyp=sqrt(dirx**2+diry**2+dirz**2)
60   - dirx=dirx/hyp
61   - diry=diry/hyp
62   - dirz=dirz/hyp
63   - costheta=dirx*posx+diry*posy+dirz*posz
64   - cond=(costheta<=1)&(costheta>=-1)&(energy>Eliminf)&(energy<Elimsup)
65   - theta = arccos(costheta[cond])*degre
66   - weight= weight[cond]
67   - energy= energy[cond]
68   - print fileId, ":", shape(energy)
69   -
70   - # Angle versus energy
71   - #=====================
72   - nbBins = 40
73   - Edata =log10(energy)
74   - angle,ener =histogram(Edata,nbBins,weights=theta)
75   - nb,ener =histogram(Edata,nbBins)
76   - ener = 10**ener
77   - enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2
78   - angle = angle/nb
79   - if argv[1] == "EGMF":
80   - ax11.plot(enercenter,angle,'.'+color[ind],
81   - label="$10^{-%.0f"%float(fileId)+"}$Gauss - MC")
82   - #ax11.plot(energy,theta,'.'+color[ind],label="$10^{-%.0f"%float(fileId)+"}$Gauss - MC")
83   - B=10**(-float(fileId))
84   - yfit = fit_theta(enercenter,B)
85   - ax11.plot(enercenter,yfit,'--'+color[ind],linewidth=2,label="$10^{-%.0f"%float(fileId)+"}$Gauss - analytic")
86   - elif argv[1] == "EBL":
87   - ax11.plot(enercenter,angle,'.'+color[ind],label=Model[int(fileId)-1])
88   - if ind==0:
89   - yfit = fit_theta(enercenter,B)
90   - ax11.plot(enercenter,yfit,'--'+color[ind],linewidth=2,label="$10^{-15}$Gauss - analytic")
91   -
92   - error=angle/yfit-1
93   - ax12.plot(enercenter,error,'+'+color[ind])
94   -
95   - # figure: radial distribution
96   - #=============================
97   - nbBins = 20
98   - #theta = theta[theta!=0]
99   - #weight=weight[theta!=0]
100   - #theta = log10(theta**2)
101   - dn,angle1=histogram(theta**2,nbBins, range=[0,0.25],weights=weight)
102   - #angle1=10**angle1
103   - thetacenter=(angle1[1:nbBins+1]+angle1[0:nbBins])/2
104   - dtheta=angle1[1:nbBins+1]-angle1[0:nbBins]
105   - dndtheta2=dn/dtheta
106   - ax2.plot(thetacenter,dndtheta2,'-'+color[ind],
107   - drawstyle='steps-mid',label="$10^{-%.0f"%float(fileId)+"}$Gauss - MC")
108   - #ax2.hist(theta, nbBins, range=[0,0.5], weights=weight,log=True,
109   - # facecolor=color[ind],alpha=.5,label="$10^{-%.0f"%float(fileId)+"}$Gauss")
110   -
111   - ind=ind+1
112   -
113   -ax11.legend(loc="best")
114   -ax11.set_xscale('log')
115   -ax11.set_yscale('log')
116   -ax11.set_ylim([0,180])
117   -ax11.grid(b=True,which='major')
118   -ax11.set_ylabel("$\\theta$ [degre]")
119   -ax12.set_xscale('log')
120   -ax12.grid(b=True,which='major')
121   -ax12.set_xlabel("energy [GeV]")
122   -ax12.set_ylabel("relative error")
123   -step=0.2*(max(error)-min(error))
124   -plt.ylim(min(error)-step, max(error)+step)
125   -xticklabels = ax11.get_xticklabels()
126   -plt.setp(xticklabels, visible=False)
127   -
128   -ax2.legend(loc="best",title="%3d - %3d GeV"%(Eliminf,Elimsup))
129   -#ax2.set_xscale('log')
130   -ax2.set_yscale('log')
131   -ax2.grid(b=True,which='major')
132   -ax2.set_xlabel("$\\theta^2$ [deg$^2$]")
133   -ax2.set_ylabel("$dN/d\\theta^2$ [deg$^{-2}$]")
134   -
135   -plt.show()
Constantes.pyc deleted
View file @e4ba931
No preview for this file type
Generation.py deleted
View file @e4ba931
... ... @@ -1,34 +0,0 @@
1   -#!/usr/bin/python
2   -
3   -from sys import argv
4   -import numpy as np
5   -import matplotlib.pyplot as plt
6   -
7   -if len(argv)<2:
8   - print "Please give at least 1 file id \n For example: 14 for Results.EGMF14 \n or 14 15 for EGMF14 and EGMF15"
9   - exit()
10   -
11   -# figure: energy distribution
12   -#=============================
13   -nbBins = 7
14   -color=['b','g','r','c','m','y']
15   -
16   -fig = plt.figure()
17   -ax = fig.add_subplot(111)
18   -
19   -ind = 0
20   -for fileId in argv[1:]:
21   - weight,nbGen=np.loadtxt("Results_EGMF"+fileId+"/Results_extra",unpack=True,usecols=[2,3])
22   -
23   - plt.hist(nbGen,nbBins,range=[2,9],log=1,facecolor=color[ind],alpha=.5, weights=weight,
24   - label="$10^{-%.0f"%float(fileId)+"}$Gauss")
25   - ind=ind+1
26   -
27   -xtick=np.arange(9)
28   -ax.legend(loc='best')
29   -ax.set_xticks(xtick+0.5)
30   -ax.set_xticklabels(xtick)
31   -ax.set_xlabel("Generation")
32   -ax.set_ylabel("Number of events")
33   -
34   -plt.show()
Map_arrival.py deleted
View file @e4ba931
... ... @@ -1,83 +0,0 @@
1   -#!/usr/bin/python
2   -
3   -from sys import argv
4   -from numpy import *
5   -import matplotlib.pyplot as plt
6   -from matplotlib.colors import LogNorm
7   -from mpl_toolkits.mplot3d.axes3d import Axes3D
8   -
9   -if len(argv)<2:
10   - print "Please give at least 1 file id \n For example: 14 for Results.EGMF14 \n or 14 15 for EGMF14 and EGMF15"
11   - exit()
12   -
13   -# figure: position of the arriving particles in space -> size circle = energy
14   -#=============================================================================
15   -nbBins = 100
16   -degre= 180/pi
17   -borne=5 # degree
18   -Elim=1e-1 # GeV
19   -
20   -fig = plt.figure()
21   -
22   -nbplots = len(argv)-1
23   -ind = 1
24   -for fileId in argv[1:]:
25   - energy,thetaDir,phiDir = loadtxt("Results_EGMF"+fileId+"/Results_momentum",unpack=True,usecols=[0,4,5])
26   - thetaPos,phiPos = loadtxt("Results_EGMF"+fileId+"/Results_position",unpack=True,usecols=[4,5])
27   - charge,weight= loadtxt("Results_EGMF"+fileId+"/Results_extra",unpack=True,usecols=[0,2])
28   -
29   - thetaPos = thetaPos*degre
30   - phiPos = phiPos*degre
31   - thetaDir = thetaDir*degre
32   - phiDir = phiDir*degre
33   - theta = thetaDir - thetaPos
34   - phi = phiDir -phiPos
35   -
36   - cond=(energy>Elim) & (abs(theta)<borne) & (abs(phi)<borne)
37   - Energy = energy[cond]
38   - weight=weight[cond]
39   - theta=theta[cond]
40   - phi=phi[cond]
41   -
42   - print "Cut at", Elim, "GeV:", fileId, "->", shape(Energy)[0], "particles"
43   - plotid = 100+nbplots*10+ind
44   -
45   - #Emax = max(Energy)
46   - #Emin = min(Energy)
47   - #area = pi * (15 * Energy/Emax)**2 +15
48   - #colors = Energy/Emax
49   -
50   - #ax = fig.add_subplot(plotid)
51   - #sax = ax.scatter(theta, phi, s=area, c=colors, alpha=0.3, marker='o')
52   - #ax.set_xlim((-borne,borne))
53   - #ax.set_ylim((-borne,borne))
54   - #ax.set_xlabel("$\\theta$ in degree")
55   - #ax.set_ylabel("$\\phi$ in degree")
56   - #ax.set_title("EGMF: $10^{-%.0f"%float(fileId)+"}$Gauss")
57   -
58   - #if Emin < 1e-3:
59   - # liminf = "%.0f"%(Emin*1e6)+"keV"
60   - #elif Emin < 1:
61   - # liminf = "%.0f"%(Emin*1e3)+"MeV"
62   - #else:
63   - # liminf = "%.0f"%(Emin)+"GeV"
64   - #limsup = "%.0f"%(Emax*1e-3)+"TeV"
65   - #cbar = fig.colorbar(sax, ticks=[0.01,1])
66   - #cbar.ax.set_yticklabels([liminf,limsup])
67   -
68   - ###########################################
69   -
70   - ax = fig.add_subplot(plotid)
71   - H,xedges,yedges,im = ax.hist2d(theta,phi,nbBins,weights=weight/shape(weight),norm=LogNorm())
72   - extent = [yedges[0], yedges[-1], xedges[0], xedges[-1]]
73   - cbar=fig.colorbar(im, ax=ax)
74   - cbar.ax.set_ylabel("counts")
75   - ax.set_xlim((-borne,borne))
76   - ax.set_ylim((-borne,borne))
77   - ax.set_xlabel("$\\theta$ in degree")
78   - ax.set_ylabel("$\\phi$ in degree")
79   - ax.set_title("EGMF: $10^{-%.0f"%float(fileId)+"}$Gauss")
80   -
81   - ind=ind+1
82   -
83   -plt.show()
Modules/Analytic.py 0 → 100644
  Diff comments   View file @7030f15
... ... @@ -0,0 +1,15 @@
  1 +from Constants import me, Ecmb, Eic, lambda_gg, degre, RL, Dic
  2 +from numpy import sqrt, abs, sin, arcsin
  3 +
  4 +def Analytic_flux(E_ic):
  5 + return me*1e3/4*sqrt(3e9/Ecmb)*1e-9*sqrt(E_ic) #GeV
  6 +
  7 +def Analytic_theta(E,fileId):
  8 + Esource=100 #TeV
  9 + Dsource=135 #Mpc
  10 + B=1e-15 #Gauss
  11 + delta_to_theta=(lambda_gg/Esource*1e3)/Dsource
  12 + RL0=RL*(Esource/2)*(1e-17/B)
  13 + Dic0=Dic/(Esource/2)
  14 + delta=Dic0/(2*RL0)*((Esource/2)**2 *Eic/E-1)
  15 + return abs(arcsin(delta_to_theta*sin(delta))*degre)
... ...
Modules/Analytic.pyc 0 → 100644
View file @7030f15
No preview for this file type
Modules/Angle.py 0 → 100644
  Diff comments   View file @7030f15
... ... @@ -0,0 +1,103 @@
  1 +from numpy import sqrt, arccos, shape, log10, histogram
  2 +from matplotlib.pyplot import figure, show, ylim, setp
  3 +from matplotlib import gridspec
  4 +from Read import ReadMomentum, ReadPosition, ReadEnergy, ReadWeight
  5 +from Analytic import Analytic_theta
  6 +from Constants import degre
  7 +
  8 +Eliminf = 1e0 #GeV
  9 +Elimsup = 1e5 #GeV
  10 +
  11 +def compute_theta(fileName):
  12 + position = ReadPosition(fileName)
  13 + momentum = ReadMomentum(fileName)
  14 + energy = ReadEnergy(fileName)
  15 + weight = ReadWeight(fileName)
  16 +
  17 + hyp=sqrt(position[0]**2+position[1]**2+position[2]**2)
  18 + position /= hyp
  19 + hyp=sqrt(momentum[0]**2+momentum[1]**2+momentum[2]**2)
  20 + momentum /= hyp
  21 +
  22 + costheta=position[0]*momentum[0]+position[1]*momentum[1]+position[2]*momentum[2]
  23 + cond=(costheta<=1)&(costheta>=-1)&(energy>Eliminf)&(energy<Elimsup)
  24 + theta = arccos(costheta[cond])*degre
  25 + weight= weight[cond]
  26 + energy= energy[cond]
  27 + print "file", fileName, "->", shape(energy)[0], "events"
  28 +
  29 + return energy, weight, theta
  30 +
  31 +def angle_vs_energy(fileName):
  32 + energy,weight,theta = compute_theta(fileName)
  33 + nbBins = 40
  34 +
  35 + energy =log10(energy)
  36 + angle,ener =histogram(energy,nbBins,weights=theta)
  37 + nb,ener =histogram(energy,nbBins)
  38 + ener = 10**ener
  39 + enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2
  40 + angle = angle/nb
  41 +
  42 + return enercenter, angle
  43 +
  44 +def drawAngle(files):
  45 + fig = figure()
  46 + gs = gridspec.GridSpec(2, 1, height_ratios=[4,1])
  47 + fig.subplots_adjust(hspace=0.001)
  48 + ax1 = fig.add_subplot(gs[0])
  49 + ax2 = fig.add_subplot(gs[1],sharex=ax1)
  50 +
  51 + for fileId in files:
  52 + energy,angle = angle_vs_energy(fileId)
  53 + p=ax1.plot(energy,angle,'.',label=fileId)
  54 +
  55 + yfit = Analytic_theta(energy,fileId)
  56 + ax1.plot(energy,yfit,'--',color=p[0].get_color(),linewidth=2,label="analytic")
  57 +
  58 + error=angle/yfit-1
  59 + ax2.plot(energy,error,'+',color=p[0].get_color())
  60 +
  61 + ax1.legend(loc="best")
  62 + ax1.set_xscale('log')
  63 + ax1.set_yscale('log')
  64 + ax1.set_ylim([0,180])
  65 + ax1.grid(b=True,which='major')
  66 + ax1.set_ylabel("$\\theta$ [degre]")
  67 + ax2.set_xscale('log')
  68 + ax2.grid(b=True,which='major')
  69 + ax2.set_xlabel("energy [GeV]")
  70 + ax2.set_ylabel("relative error")
  71 + step=0.2*(max(error)-min(error))
  72 + ylim(min(error)-step, max(error)+step)
  73 + xticklabels = ax1.get_xticklabels()
  74 + setp(xticklabels, visible=False)
  75 +
  76 + show()
  77 +
  78 +def radial(fileName):
  79 + nbBins = 20
  80 + energy,weight,theta = compute_theta(fileName)
  81 +
  82 + dn,angle=histogram(theta**2,nbBins,range=[0,0.25],weights=weight)
  83 + thetacenter=(angle[1:nbBins+1]+angle[0:nbBins])/2
  84 + dtheta=angle[1:nbBins+1]-angle[0:nbBins]
  85 + dndtheta2=dn/dtheta
  86 +
  87 + return thetacenter, dndtheta2
  88 +
  89 +def drawRadial(files):
  90 + fig = figure()
  91 + ax = fig.add_subplot(111)
  92 +
  93 + for fileId in files:
  94 + theta,dndtheta2=radial(fileId)
  95 + ax.plot(theta,dndtheta2,drawstyle='steps-mid',label=fileId)
  96 +
  97 + ax.legend(loc="best",title="%3d - %3d GeV"%(Eliminf,Elimsup))
  98 + ax.set_yscale('log')
  99 + ax.grid(b=True,which='major')
  100 + ax.set_xlabel("$\\theta^2$ [deg$^2$]")
  101 + ax.set_ylabel("$dN/d\\theta^2$ [deg$^{-2}$]")
  102 +
  103 + show()
... ...
Modules/Angle.pyc 0 → 100644
View file @7030f15
No preview for this file type
Constantes.py renamed to Modules/Constants.py
  Diff comments   View file @7030f15
1 1 #!/bin/python
2 2  
3   -from numpy import *
4   -from scipy.special import *
  3 +from numpy import pi, sqrt
  4 +from scipy.special import zeta
5 5  
6 6 e=4.8032068e-10 # esu (cgs units)
7 7 qe=1.602176565e-12 # erg
... ... @@ -67,13 +67,3 @@ Ethreshold_ic = Eic*Dic/RL # GeV
67 67 Ethreshold_gg=(me)**2/Eebl *1e-3 #GeV
68 68 lambda_gg=0.8 # (E_gamma/1TeV)^-1 Gpc (from Durrer and Neronov 2013)
69 69  
70   -# usefull integrand
71   -def comobileTime(z):
72   - return -1/(H0*(1+z)*sqrt(omegaM*(1+z)**3+omegaK*(1+z)**2+omegaL))
73   -
74   -def distPhoton(z):
75   - return -c/(H0*a0*sqrt(omegaM*(1+z)**3+omegaK*(1+z)**2+omegaL))
76   -
77   -def distLepton(z,E):
78   - beta = sqrt(1-m**2*c**4/E**2)
79   - return -beta*c/(H0*a0*sqrt(omegaM*(1+z)**3+omegaK*(1+z)**2+omegaL))
... ...
Modules/Constants.pyc 0 → 100644
View file @7030f15
No preview for this file type
Modules/Generation.py 0 → 100644
  Diff comments   View file @7030f15
... ... @@ -0,0 +1,23 @@
  1 +from numpy import shape, arange
  2 +from matplotlib.pyplot import figure, show, hist
  3 +from Read import ReadGeneration, ReadWeight
  4 +
  5 +
  6 +def drawGeneration(files):
  7 + fig = figure()
  8 + ax = fig.add_subplot(111)
  9 + nbBins = 7
  10 + for fileId in files:
  11 + weight = ReadWeight(fileId)
  12 + generation = ReadGeneration(fileId)
  13 + print "file", fileId, "->", shape(generation)[0], "events"
  14 + hist(generation,nbBins,range=[2,9],log=1,alpha=.5, weights=weight,label=fileId)
  15 +
  16 + xtick=arange(9)
  17 + ax.legend(loc='best')
  18 + ax.set_xticks(xtick+0.5)
  19 + ax.set_xticklabels(xtick)
  20 + ax.set_xlabel("Generation")
  21 + ax.set_ylabel("Number of events")
  22 +
  23 + show()
... ...
Modules/Generation.pyc 0 → 100644
View file @7030f15
No preview for this file type
Modules/Integrand.py 0 → 100644
  Diff comments   View file @7030f15
... ... @@ -0,0 +1,12 @@
  1 +from Constants import H0, omegaK, omegaM, omegaL, c, m , a0
  2 +from numpy import sqrt
  3 +
  4 +def comobileTime(z):
  5 + return -1/(H0*(1+z)*sqrt(omegaM*(1+z)**3+omegaK*(1+z)**2+omegaL))
  6 +
  7 +def distPhoton(z):
  8 + return -c/(H0*a0*sqrt(omegaM*(1+z)**3+omegaK*(1+z)**2+omegaL))
  9 +
  10 +def distLepton(z,E):
  11 + beta = sqrt(1-m**2*c**4/E**2)
  12 + return -beta*c/(H0*a0*sqrt(omegaM*(1+z)**3+omegaK*(1+z)**2+omegaL))
... ...
Modules/Integrand.pyc 0 → 100644
View file @7030f15
No preview for this file type
Modules/Map.py 0 → 100644
  Diff comments   View file @7030f15
... ... @@ -0,0 +1,88 @@
  1 +from numpy import shape, pi
  2 +from matplotlib.pyplot import figure, show
  3 +from matplotlib.colors import LogNorm
  4 +from Read import ReadEnergy, ReadWeight, ReadMomentumAngle, ReadPositionAngle
  5 +from Constants import degre
  6 +
  7 +Elim=1e-1 # GeV
  8 +borne=5 # degree
  9 +nbBins = 100
  10 +
  11 +def thetaphi(fileId):
  12 + energy = ReadEnergy(fileId)
  13 + weight = ReadWeight(fileId)
  14 + thetaDir,phiDir = ReadMomentumAngle(fileId)*degre
  15 + thetaPos,phiPos = ReadPositionAngle(fileId)*degre
  16 +
  17 + theta = thetaDir - thetaPos
  18 + phi = phiDir -phiPos
  19 +
  20 + cond=(energy>Elim) & (abs(theta)<borne) & (abs(phi)<borne)
  21 + energy = energy[cond]
  22 + weight=weight[cond]
  23 + theta=theta[cond]
  24 + phi=phi[cond]
  25 +
  26 + print "file", fileId, "->", shape(energy)[0], "events"
  27 + return energy, theta, phi, weight
  28 +
  29 +def drawMap(files):
  30 + print "Cut at", Elim, "GeV"
  31 + fig = figure()
  32 + nbplots = len(files)
  33 +
  34 + ind=1
  35 + for fileId in files:
  36 + energy, theta, phi, weight = thetaphi(fileId)
  37 +
  38 + ax = fig.add_subplot(100+nbplots*10+ind)
  39 +
  40 + Emax = max(energy)
  41 + Emin = min(energy)
  42 + colors = energy/Emax
  43 + area = pi * (15 * colors)**2 +15
  44 + sax = ax.scatter(theta, phi, c=colors, s=area, alpha=0.3, marker='o')
  45 +
  46 + if Emin < 1e-3:
  47 + liminf = "%.0f"%(Emin*1e6)+"keV"
  48 + elif Emin < 1:
  49 + liminf = "%.0f"%(Emin*1e3)+"MeV"
  50 + else:
  51 + liminf = "%.0f"%(Emin)+"GeV"
  52 + limsup = "%.0f"%(Emax*1e-3)+"TeV"
  53 + cbar = fig.colorbar(sax, ticks=[0.01,1])
  54 + cbar.ax.set_yticklabels([liminf,limsup])
  55 +
  56 + ax.set_xlim((-borne,borne))
  57 + ax.set_ylim((-borne,borne))
  58 + ax.set_xlabel("$\\theta$ [deg]")
  59 + ax.set_ylabel("$\\phi$ [deg]")
  60 + ax.set_title(fileId)
  61 +
  62 + ind=ind+1
  63 +
  64 + show()
  65 +
  66 +def drawMapHisto(files):
  67 + print "Cut at", Elim, "GeV"
  68 + fig = figure()
  69 + nbplots = len(files)
  70 +
  71 + ind=1
  72 + for fileId in files:
  73 + energy, theta, phi, weight = thetaphi(fileId)
  74 +
  75 + ax = fig.add_subplot(100+nbplots*10+ind)
  76 + H,xedges,yedges,im = ax.hist2d(theta,phi,nbBins,weights=weight/shape(weight),norm=LogNorm())
  77 + extent = [yedges[0], yedges[-1], xedges[0], xedges[-1]]
  78 + cbar=fig.colorbar(im, ax=ax)
  79 + cbar.ax.set_ylabel("counts")
  80 + ax.set_xlim((-borne,borne))
  81 + ax.set_ylim((-borne,borne))
  82 + ax.set_xlabel("$\\theta$ [deg]")
  83 + ax.set_ylabel("$\\phi$ [deg]")
  84 + ax.set_title(fileId)
  85 +
  86 + ind=ind+1
  87 +
  88 + show()
... ...
Modules/Map.pyc 0 → 100644
View file @7030f15
No preview for this file type
Modules/Read.py 0 → 100644
  Diff comments   View file @7030f15
... ... @@ -0,0 +1,44 @@
  1 +from numpy import loadtxt
  2 +
  3 +resultDirectory="Results/"
  4 +positionFile="/position"
  5 +momentumFile="/momentum"
  6 +extraFile="/extra"
  7 +ElmagFile="/spec_diff_test"
  8 +
  9 +def ReadTime(fileName):
  10 + return loadtxt(resultDirectory+fileName+positionFile,unpack=True,usecols=[0])
  11 +
  12 +def ReadPosition(fileName):
  13 + return loadtxt(resultDirectory+fileName+positionFile,unpack=True,usecols=[1,2,3])
  14 +
  15 +def ReadPositionAngle(fileName):
  16 + return loadtxt(resultDirectory+fileName+positionFile,unpack=True,usecols=[4,5])
  17 +
  18 +def ReadEnergy(fileName):
  19 + return loadtxt(resultDirectory+fileName+momentumFile,unpack=True,usecols=[0])
  20 +
  21 +def ReadMomentum(fileName):
  22 + return loadtxt(resultDirectory+fileName+momentumFile,unpack=True,usecols=[1,2,3])
  23 +
  24 +def ReadMomentumAngle(fileName):
  25 + return loadtxt(resultDirectory+fileName+momentumFile,unpack=True,usecols=[4,5])
  26 +
  27 +def ReadCharge(fileName):
  28 + return loadtxt(resultDirectory+fileName+extraFile,unpack=True,usecols=[0])
  29 +
  30 +def ReadRedshift(fileName):
  31 + return loadtxt(resultDirectory+fileName+extraFile,unpack=True,usecols=[1])
  32 +
  33 +def ReadWeight(fileName):
  34 + return loadtxt(resultDirectory+fileName+extraFile,unpack=True,usecols=[2])
  35 +
  36 +def ReadGeneration(fileName):
  37 + return loadtxt(resultDirectory+fileName+extraFile,unpack=True,usecols=[3])
  38 +
  39 +def ReadElmag(fileName):
  40 + Elmag=resultDirectory+fileName+ElmagFile
  41 + energy,fluxGamma,fluxLepton = loadtxt(Elmag,unpack=True,usecols=[0,1,2])
  42 + energy *= 10**(-9)
  43 + flux=fluxGamma+fluxLepton
  44 + return energy,flux
... ...
Modules/Read.pyc 0 → 100644
View file @7030f15
No preview for this file type
Modules/Spectrum.py 0 → 100644
  Diff comments   View file @7030f15
... ... @@ -0,0 +1,45 @@
  1 +from numpy import histogram, log10, shape, logspace
  2 +from matplotlib.pyplot import figure, show
  3 +from os.path import isfile
  4 +from Read import ReadEnergy, ReadWeight, ReadElmag, resultDirectory, ElmagFile
  5 +from Analytic import Analytic_flux
  6 +
  7 +def spectrum(fileName):
  8 + energy = ReadEnergy(fileName)
  9 + weight = ReadWeight(fileName)
  10 + print "file", fileName, "->", shape(energy)[0], "events"
  11 + nbBins=100
  12 + energy =log10(energy)
  13 + flux,ener=histogram(energy,nbBins,weights=weight)
  14 + ener = 10**ener
  15 + enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2
  16 + binSize=ener[1:nbBins+1]-ener[0:nbBins]
  17 + flux = enercenter**2 * flux/binSize
  18 + return enercenter,flux
  19 +
  20 +def drawSpectrum(files):
  21 + fig = figure()
  22 + ax = fig.add_subplot(111)
  23 + for fileId in files:
  24 + energy,flux = spectrum(fileId)
  25 + maxflux = 10 #shape(energy)[0]
  26 + flux /= maxflux
  27 + p=ax.plot(energy,flux,drawstyle='steps-mid',label=fileId)
  28 +
  29 + Elmag=resultDirectory+fileId+ElmagFile
  30 + if isfile(Elmag):
  31 + energy, flux = ReadElmag(fileId)
  32 + ax.plot(energy,flux,"-",color=p[0].get_color(),label="Elmag - "+fileId)
  33 +
  34 + xfit = logspace(log10(10**(-3)),log10(2*10**4),200)
  35 + yfit = Analytic_flux(xfit)
  36 + ax.plot(xfit,yfit,'--m',linewidth=2,label="Analytic ($\\sqrt{E}$)")
  37 +
  38 + ax.set_xscale('log')
  39 + ax.set_yscale('log')
  40 + ax.grid(b=True,which='major')
  41 + ax.legend(loc="best")
  42 + ax.set_xlabel("energy [GeV]")
  43 + ax.set_ylabel("$E^2 \\frac{dN}{dE}$ [GeV]")
  44 +
  45 + show()
... ...
Modules/Spectrum.pyc 0 → 100644
View file @7030f15
No preview for this file type
Modules/Timing.py 0 → 100644
  Diff comments   View file @7030f15
... ... @@ -0,0 +1,43 @@
  1 +from numpy import histogram, log10, shape
  2 +from matplotlib.pyplot import figure, show
  3 +from scipy.integrate import quad
  4 +from Read import ReadTime, ReadWeight
  5 +from Constants import yr
  6 +from Integrand import comobileTime
  7 +
  8 +def timing(fileName):
  9 + time=ReadTime(fileName)
  10 + weight=ReadWeight(fileName)
  11 + prop = float(shape(time[time<0])[0])/shape(time)[0]
  12 + print "file", fileName, "->", shape(time)[0], "events", "negative time:", shape(time[time<0]), "~", prop
  13 +
  14 + time=time[time>0]
  15 + weight=weight[time>0]
  16 +
  17 + nbBins = 100
  18 + zSource = 0.0308
  19 + tlim = quad(comobileTime,zSource,0)[0]/yr
  20 + time=log10(time/tlim)
  21 +
  22 + dN,dt=histogram(time,nbBins,weights=weight)
  23 + timecenter=(dt[1:nbBins+1]+dt[0:nbBins])/2
  24 + binSize=dt[1:nbBins+1]-dt[0:nbBins]
  25 + Nmax=shape(time)[0]
  26 + dNdt=dN/(Nmax*binSize)
  27 +
  28 + return timecenter, dNdt
  29 +
  30 +def drawTiming(files):
  31 + fig = figure()
  32 + ax = fig.add_subplot(111)
  33 + for fileId in files:
  34 + time,dNdt = timing(fileId)
  35 + ax.plot(time[dNdt!=0],dNdt[dNdt!=0],linestyle="steps-mid",label=fileId)
  36 +
  37 + ax.set_yscale('log')
  38 + ax.grid(b=True,which='major')
  39 + ax.legend(loc="best")
  40 + ax.set_xlabel("Time [$\Delta t / t$ log scale]")
  41 + ax.set_ylabel("$t\ dN/dt$")
  42 +
  43 + show()
... ...
Modules/Timing.pyc 0 → 100644
View file @7030f15
No preview for this file type
Modules/__init__.py 0 → 100644
  Diff comments   View file @7030f15
Modules/__init__.pyc 0 → 100644
View file @7030f15
No preview for this file type
Spectrum.py deleted
View file @e4ba931
... ... @@ -1,87 +0,0 @@
1   -#!/usr/bin/python
2   -
3   -from sys import argv
4   -import matplotlib.pyplot as plt
5   -from Constantes import *
6   -
7   -if argv[1] == "EGMF":
8   - fileName="Results_EGMF"
9   -elif argv[1] == "EBL":
10   - fileName="Results_EBL"
11   - B=10**(-15)
12   - Model=["Kneiske and Doll - 'best fit'","Kneiske and Doll - 'lower limit'",
13   - "Franceschini","Finke and Al","Gilmore and Al","Dominguez and Al"]
14   -else:
15   - print "Used: ./Energy_distribution.py arg1 ind1 ind2 ..."
16   - print " arg1 = EGMF or EBL "
17   - print " ind1 = file number \n"
18   - print "Examples: "
19   - print " to study EGMF14: ./Energy_distribution.py EGMF 14 "
20   - print " to study EBL1 and EBL2: ./Energy_distribution.py EBL 1 2 \n"
21   - exit()
22   -
23   -# figure: energy distribution
24   -#=============================
25   -nbBins = 100
26   -convert= 180/pi
27   -color=['b','r','g','c','m','y']
28   -
29   -fig = plt.figure()
30   -ax = fig.add_subplot(111)
31   -
32   -energy= loadtxt(fileName+argv[2]+"/Results_momentum",unpack=True,usecols=[0])
33   -Edata =log10(energy)
34   -Emax = max(Edata)
35   -Emin = min(Edata)
36   -
37   -ind = 0
38   -for fileId in argv[2:]:
39   - energy= loadtxt(fileName+fileId+"/Results_momentum",unpack=True,usecols=[0])
40   - charge,weight=loadtxt(fileName+fileId+"/Results_extra",unpack=True,usecols=[0,2])
41   -
42   - cond= (charge==0)
43   - energy=energy[cond]
44   - weight=weight[cond]
45   - print "data shape:", shape(energy)
46   -
47   - Edata =log10(energy)
48   - flux,ener=histogram(Edata,nbBins,range=(Emin,Emax),weights=weight)
49   - ener = 10**ener
50   - enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2
51   - binSize=ener[1:nbBins+1]-ener[0:nbBins]
52   - flux = enercenter**2 * flux/binSize
53   -
54   - maxflux = shape(energy)[0]
55   -
56   - if argv[1] == "EGMF":
57   - ax.plot(enercenter,flux[:]/maxflux,color=color[ind],drawstyle='steps-mid',
58   - label="$J_\gamma$ $10^{-%.0f"%float(fileId)+"}$Gauss")
59   - elif argv[1] == "EBL":
60   - ax.plot(enercenter,flux[:]/maxflux,color=color[ind],drawstyle='steps-mid',
61   - label=Model[int(fileId)-1])
62   - # Results from Elmag
63   - Elmag=fileName+fileId+"/spec_diff_test"
64   - if os.path.isfile(Elmag):
65   - energy,fluxGamma,fluxLepton = loadtxt(Elmag,unpack=True,usecols=[0,1,2])
66   - energy=energy*10**(-9)
67   - flux=fluxGamma+fluxLepton
68   - flux=flux*2.5e-12
69   - ax.plot(energy,flux,"-"+color[ind],label="Elmag - "+Model[int(fileId)-1])
70   -
71   - ind=ind+1
72   -
73   -def fit(E_ic):
74   - return me*1e3/4*sqrt(3e9/Ecmb)*1e-9*sqrt(E_ic) #GeV
75   -
76   -xfit = logspace(log10(10**(-3)),log10(2*10**4),200)
77   -yfit = fit(xfit)
78   -ax.plot(xfit,yfit,'--m',linewidth=2,label="$E^{0.5}$")
79   -
80   -ax.set_xscale('log')
81   -ax.set_yscale('log')
82   -ax.grid(b=True,which='major')
83   -ax.legend(loc="best")
84   -ax.set_xlabel("energy [GeV]")
85   -ax.set_ylabel("$E^2 \\frac{dN}{dE}$ [GeV.photon]")
86   -
87   -plt.show()
Time_delay.py deleted
View file @e4ba931
... ... @@ -1,96 +0,0 @@
1   -#!/usr/bin/python
2   -
3   -from sys import argv
4   -import os.path
5   -import matplotlib.pyplot as plt
6   -from Constantes import *
7   -from scipy.integrate import quad
8   -
9   -if argv[1] == "EGMF":
10   - fileName="Results_EGMF"
11   -elif argv[1] == "EBL":
12   - fileName="Results_EBL"
13   - B=10**(-15)
14   - Model=["Kneiske and Doll - 'best fit'","Kneiske and Doll - 'lower limit'",
15   - "Franceschini","Finke and Al","Gilmore and Al","Dominguez and Al"]
16   -elif argv[1] == "prec":
17   - fileName="Results_prec"
18   - B=10**(-17)
19   -else:
20   - print "Used: ./Energy_distribution.py arg1 ind1 ind2 ..."
21   - print " arg1 = EGMF or EBL "
22   - print " ind1 = file number \n"
23   - print "Examples: "
24   - print " to study EGMF14: ./Energy_distribution.py EGMF 14 "
25   - print " to study EBL1 and EBL2: ./Energy_distribution.py EBL 1 2 \n"
26   - exit()
27   -
28   -# figure: energy distribution
29   -#=============================
30   -zSource = 0.0308
31   -nbBins = 100
32   -convert= 180/pi
33   -color=['b','r','g','c','m','y']
34   -
35   -fig = plt.figure()
36   -ax = fig.add_subplot(111)
37   -
38   -tlim = quad(comobileTime,zSource,0)[0]/yr
39   -print tlim
40   -
41   -time= loadtxt(fileName+argv[2]+"/Results_position",unpack=True,usecols=[0])
42   -time=time[time>0]
43   -time=time/tlim
44   -tmax = max(time)
45   -tmin = min(time)
46   -print tmin, tmax
47   -
48   -ind = 0
49   -for fileId in argv[2:]:
50   - time= loadtxt(fileName+fileId+"/Results_position",unpack=True,usecols=[0])
51   - charge,weight,gen=loadtxt(fileName+fileId+"/Results_extra",unpack=True,usecols=[0,2,3])
52   - time=time/tlim
53   -
54   - prop = float(shape(time[time<0])[0])/shape(time)[0]
55   - print "data shape",fileName+fileId,":", shape(time), "negative time:", shape(time[time<0]), "~", prop
56   -
57   - cond= (gen<10) & (time>0) & (charge==0)
58   - time=time[cond]
59   - weight=weight[cond]
60   - time=log10(time)
61   -
62   - dN,dt=histogram(time,nbBins,range=(log10(tmin),log10(tmax)),weights=weight)
63   - #time=10**time
64   - timecenter=(dt[1:nbBins+1]+dt[0:nbBins])/2
65   - binSize=dt[1:nbBins+1]-dt[0:nbBins]
66   - dNdt=dN/binSize
67   -
68   - Nmax=shape(time)[0]
69   -
70   - if argv[1] == "EGMF":
71   - #ax.hist(time,nbBins,weights=weight,range=(log10(tmin),log10(tmax)),log=1,facecolor=color[ind],
72   - # alpha=.5, label="$10^{-%.0f"%float(fileId)+"}$Gauss")
73   - ax.plot(timecenter[dNdt!=0],dNdt[dNdt!=0]/Nmax,color=color[ind],linestyle="steps-mid",
74   - label="$10^{-%.0f"%float(fileId)+"}$Gauss")
75   - elif argv[1] == "prec":
76   - ax.hist(time,nbBins,weights=weight,range=(log10(tmin),log10(tmax)),log=1,facecolor=color[ind],
77   - alpha=.5)#, label="prec = $10^{-%.0f"%float(fileId)+"}$")
78   - ax.plot(timecenter[dNdt!=0],dNdt[dNdt!=0],"-"+color[ind],
79   - label="prec = $10^{-%.0f"%float(fileId)+"}$")
80   - elif argv[1] == "EBL":
81   - #plt.hist(time,nbBins,weights=weight,range=(tmin,tmax),log=1,facecolor=color[ind],alpha=.5,
82   - ax.plot(timecenter,dNdt,"."+color[ind],
83   - label=Model[int(fileId)-1])
84   -
85   - ind=ind+1
86   -
87   -#ax.set_xscale('log')
88   -ax.set_yscale('log')
89   -ax.grid(b=True,which='major')
90   -ax.legend(loc="best")
91   -if argv[1]=="prec":
92   - ax.legend(loc="best",title="$10^{-17}$Gauss")
93   -ax.set_xlabel("Time [$\Delta t / t$ log scale]")
94   -ax.set_ylabel("$t\ dN/dt$")
95   -
96   -plt.show()