Thomas Fitoussi / simulation-analysis

Blame view

Angle_distribution.py 4.04 KB
Edit Raw Normal View History
b5217158   Thomas Fitoussi   Python scripts to... 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
 
#!/usr/bin/python

from sys import argv
import matplotlib.pyplot as plt
from matplotlib import gridspec
from Constantes import *

if argv[1] == "EGMF":
   fileName="Results_EGMF"
elif argv[1] == "EBL":
   fileName="Results_EBL"
   B=10**(-15)
   Model=["Kneiske and Doll - 'best fit'","Kneiske and Doll - 'lower limit'",
         "Fraceschini","Finke and Al","Gilmore and Al","Dominguez and Al"]
else:
   print "Used: ./Angle_distribution.py arg1 ind1 ind2 ..."
   print "        arg1 = EGMF or EBL "
   print "        ind1 = file number \n"
   print "Examples: "
   print "   to study EGMF14: ./Angle_distribution.py EGMF 14 "
   print "   to study EBL1 and EBL2: ./Angle_distribution.py EBL 1 2 \n"
   exit()

nbBins = 100

Esource=100 #TeV
Dsource=135 #Mpc
delta_to_theta=(lambda_gg/Esource*1e3)/Dsource

def fit_theta(E,B):
   RL0=RL*(Esource/2)*(1e-17/B)
   Dic0=Dic/(Esource/2)
   delta=Dic0/(2*RL0)*((Esource/2)**2 *Eic/E-1)
   return abs(arcsin(delta_to_theta*sin(delta))*degre)

#fig2 = plt.figure()
#ax2 = fig2.add_subplot(111)

fig1 = plt.figure()
gs = gridspec.GridSpec(2, 1, height_ratios=[4,1]) 
fig1.subplots_adjust(hspace=0.001)
ax11 = fig1.add_subplot(gs[0])
ax12 = fig1.add_subplot(gs[1],sharex=ax11)

color=['b','r','g','c','m','y']

ind = 0
for fileId in argv[2:]:
   energy,dirx,diry,dirz = loadtxt(fileName+fileId+"/Results_momentum",unpack=True,usecols=[0,1,2,3])
   posx,posy,posz = loadtxt(fileName+fileId+"/Results_position",unpack=True,usecols=[1,2,3])
4448ea24   Thomas Fitoussi   Time distribution 
51
 
   charge,weight,gen=loadtxt(fileName+fileId+"/Results_extra",unpack=True,usecols=[0,2,3])
b5217158   Thomas Fitoussi   Python scripts to... 
52
53
54
55
56
57
58
59
60
61
 

   hyp=sqrt(posx**2+posy**2+posz**2)
   posx=posx/hyp
   posy=posy/hyp
   posz=posz/hyp
   hyp=sqrt(dirx**2+diry**2+dirz**2)
   dirx=dirx/hyp
   diry=diry/hyp
   dirz=dirz/hyp
   costheta=dirx*posx+diry*posy+dirz*posz
4448ea24   Thomas Fitoussi   Time distribution 
62
 
   cond=(costheta<=1)&(costheta>=-1)&(charge==0)
b5217158   Thomas Fitoussi   Python scripts to... 
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
 
   theta = arccos(costheta[cond])*degre
   weight= weight[cond]
   energy= energy[cond]
   print fileId, ":", shape(energy)

   # Angle versus energy
   #===================== 
   nbBins = 40
   Edata =log10(energy)
   angle,ener =histogram(Edata,nbBins,weights=theta)
   nb,ener =histogram(Edata,nbBins)
   ener = 10**ener
   enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2
   angle = angle/nb 
   if argv[1] == "EGMF":
      ax11.plot(enercenter,angle,'.'+color[ind],
               label="$10^{-%.0f"%float(fileId)+"}$Gauss - MC")
      #ax11.plot(energy,theta,'.'+color[ind],label="$10^{-%.0f"%float(fileId)+"}$Gauss - MC")
      B=10**(-float(fileId))
      yfit = fit_theta(enercenter,B)
      ax11.plot(enercenter,yfit,'--'+color[ind+1],linewidth=2,label="$10^{-%.0f"%float(fileId)+"}$Gauss - analytic")
   elif argv[1] == "EBL":
      ax11.plot(enercenter,angle,'.'+color[ind],label=Model[int(fileId)-1])
      if ind==0: 
         yfit = fit_theta(enercenter,B)
         ax11.plot(enercenter,yfit,'--'+color[ind+1],linewidth=2,label="$10^{-15}$Gauss - analytic")

   error=angle/yfit-1
   ax12.plot(enercenter,error,'+'+color[ind])

   # figure: radial distribution
   #=============================
   #nbBins = 20
   #theta = log10(theta)
   #nb,angle1=histogram(theta, nbBins, weights=weight)
   #angle1=10**angle1
   #anglecenter=(angle1[1:nbBins+1]+angle1[0:nbBins])/2
   #binSize=angle1[1:nbBins+1]-angle1[0:nbBins]
   #ax2.plot(anglecenter,nb,'-'+color[ind],label="$10^{-%.0f"%float(fileId)+"}$Gauss - MC")
   #ax2.hist(theta, nbBins, weights=weight,log=True,
   #      facecolor=color[ind],alpha=.5,label="$10^{-%.0f"%float(fileId)+"}$Gauss")

   ind=ind+1

ax11.legend(loc="best")
ax11.set_xscale('log')
ax11.set_yscale('log')
ax11.set_ylim([0,180])
ax11.grid(b=True,which='major')
ax11.set_ylabel("$\\theta$ [degre]")
ax12.set_xscale('log')
ax12.grid(b=True,which='major')
ax12.set_xlabel("energy [GeV]")
ax12.set_ylabel("relative error")
step=0.2*(max(error)-min(error))
plt.ylim(min(error)-step, max(error)+step)
xticklabels = ax11.get_xticklabels()
plt.setp(xticklabels, visible=False)

#ax2.legend(loc="best")
#ax2.set_xscale('log')
#ax2.set_yscale('log')
#ax2.grid(b=True,which='major')
#ax2.set_xlabel("$\\theta$ [degre]")
#ax2.set_ylabel("Nb events * weight")

plt.show()