Thomas Fitoussi / simulation-analysis

Commit 0b56ded50cdb37ff0f69daba43d96cd12425405e

Authored by Thomas Fitoussi
1 parent 28843f23
Exists in master

take into account source angle of emission 

=> selection for a cone emission rather than isotropic
Showing 5 changed files with 118 additions and 92 deletions   Show diff stats
Inline Side-by-side
Analysis.py
  Diff comments   View file @0b56ded
... ... @@ -5,7 +5,7 @@ from numpy import append, savetxt, shape, array, newaxis, zeros, arange
5 5 from Modules.Read import ReadEnergy, ReadTime, ReadExtraFile, ReadProfile
6 6 from Modules.Read import ReadMomentumAngle, ReadPositionAngle, resultDirectory
7 7 from Modules.Spectrum import spectrum
8   -from Modules.Map import computeMap, isotrop
  8 +from Modules.Map import computeMap
9 9 from Modules.Angle import angle_vs_energy
10 10 from Modules.Timing import timing
11 11 from Modules.Constants import degre
... ... @@ -19,6 +19,7 @@ if shape(argv)[0] < 2:
19 19  
20 20 #=============================================================================#
21 21 PowerSpectrum=[1,1.5,2,2.5]
  22 +Jet_Opening=[180,30,60] # degre (180 <=> isotrop)
22 23  
23 24 for fileId in argv[1:]:
24 25 print "#=============================================================================#"
... ... @@ -28,14 +29,14 @@ for fileId in argv[1:]:
28 29 print "# 1. Reading data"
29 30 time = ReadTime(fileId)
30 31 energy = ReadEnergy(fileId)
31   - weightini, generation, theta_arrival, Esource = ReadExtraFile(fileId,[2,3,4,5])
  32 + weightini, generation, theta_arrival, Esource, dir_source = ReadExtraFile(fileId,[2,3,4,5,6])
32 33 nbPhotonsEmitted=ReadProfile(fileId,[3])
  34 + weightini /= nbPhotonsEmitted
33 35  
34 36 Gen_contrib = []
35 37 Gen_cont = zeros((int(max(generation))+1))
36 38 Source = []
37 39 Spectrum = []
38   - Angle_Energy = []
39 40 Timing = []
40 41  
41 42 print "# 2. Computing powerlaw spectrum"
... ... @@ -61,7 +62,7 @@ for fileId in argv[1:]:
61 62 nbBins = 50
62 63 # draw source spectrum
63 64 Es=array(list(set(Esource)))
64   - Ws= (Es/min(Es))**(1-powerlaw_index)
  65 + Ws= (Es/min(Es))**(1-powerlaw_index) /nbPhotonsEmitted
65 66 Es,Fs = spectrum(Es,Ws,nbBins=nbBins)
66 67 Es=Es[:,newaxis]
67 68 Fs=Fs[:,newaxis]
... ... @@ -102,44 +103,42 @@ for fileId in argv[1:]:
102 103 #=============================================================================#
103 104 # IMAGING, RADIAL DISTRIBUTION AND ANGLE VERSUS ENERGY
104 105 #=============================================================================#
105   - print "# 3. compute images and radial distribution"
106   - powerlaw_index = 2
107   - Elim = 1e-1 # GeV
108   - thetalim = 20 # degre
109   - nbBins = 100
  106 + Angle_Energy = []
110 107 Radial = []
  108 + print "# 3. compute images and radial distribution"
  109 + for jet_opening_angle in Jet_Opening:
  110 + powerlaw_index = 2
  111 + Elim = 1e-1 # GeV
  112 + thetalim = 20 # degre
  113 + nbBins = 100
  114 +
  115 + thetaDir,phiDir = ReadMomentumAngle(fileId)*degre
  116 + thetaPos,phiPos = ReadPositionAngle(fileId)*degre
  117 + theta = thetaDir - thetaPos
  118 + phi = phiDir -phiPos
  119 +
  120 + # apply source spectrum
  121 + weight_source = (Esource/min(Esource))**(1-powerlaw_index)
  122 + weight = weightini* weight_source
  123 +
  124 + # apply selection
  125 + theta2,dndtheta2,ener,angle = computeMap(theta,phi,weight,energy,dir_source,fileId,
  126 + jet_opening_angle,Elim,nbBins,borne=[thetalim,thetalim])
  127 + theta2=theta2[:,newaxis]
  128 + dndtheta2=dndtheta2[:,newaxis]
  129 + ener=ener[:,newaxis]
  130 + angle=angle[:,newaxis]
  131 +
  132 + if Radial==[]:
  133 + Radial = theta2
  134 + Radial = append(Radial,dndtheta2,axis=1)
  135 + Angle_Energy = ener
  136 + Angle_Energy = append(Angle_Energy,angle,axis=1)
  137 + else:
  138 + Radial = append(Radial,dndtheta2,axis=1)
  139 + Angle_Energy = append(Angle_Energy,angle,axis=1)
111 140  
112   - thetaDir,phiDir = ReadMomentumAngle(fileId)*degre
113   - thetaPos,phiPos = ReadPositionAngle(fileId)*degre
114   - theta = thetaDir - thetaPos
115   - phi = phiDir -phiPos
116   -
117   - # apply source spectrum
118   - weight_source = (Esource/min(Esource))**(1-powerlaw_index)
119   - weight = weightini* weight_source
120   -
121   - # apply selection
122   - cond=(energy>Elim) #& (abs(theta)<thetalim) & (abs(phi)<thetalim)
123   - theta=theta[cond]
124   - phi=phi[cond]
125   - weight=weight[cond]
126   - energy=energy[cond]
127   - theta2,dndtheta2 = computeMap(theta,phi,weight,energy,fileId,nbBins,
128   - source=isotrop(),borne=[thetalim,thetalim])
129   - ener,angle = angle_vs_energy(theta,energy,weight,nbBins)
130   - theta2=theta2[:,newaxis]
131   - dndtheta2=dndtheta2[:,newaxis]
132   - ener=ener[:,newaxis]
133   - angle=angle[:,newaxis]
134   -
135   - if Radial==[]:
136   - Radial = theta2
137   - Radial = append(Radial,dndtheta2,axis=1)
138   - Angle_Energy = ener
139   - Angle_Energy = append(Angle_Energy,angle,axis=1)
140   - else:
141   - Radial = append(Radial,dndtheta2,axis=1)
142   - Angle_Energy = append(Angle_Energy,angle,axis=1)
  141 + InProgress(Jet_Opening.index(jet_opening_angle)+1,shape(Jet_Opening)[0])
143 142  
144 143 #=============================================================================#
145 144 print "# 4. writing files"
... ...
Modules/Map.py
  Diff comments   View file @0b56ded
1 1 from numpy import shape, histogram2d, arange, newaxis, array, nditer, zeros, histogram
2 2 from numpy import sqrt, pi, exp, log, log10, sin, cos, where
3   -from scipy.signal import convolve2d
4 3 from matplotlib.pyplot import figure, matshow, savefig, show
5 4 from matplotlib.colors import LogNorm
6 5 from Read import resultDirectory, ReadRadial, ReadAngleVsEnergy
... ... @@ -15,16 +14,7 @@ philim=float(philim)
15 14 thetalim=float(thetalim)
16 15 limits = [[-thetalim,thetalim],[-philim,philim]]
17 16  
18   -def isotrop():
19   - return array([[1]]), "isotropic"
20   -
21   -def Gaussian2D(center=[0,0],fwhm=0.01):
22   - x = arange(-thetalim,thetalim+step,step)
23   - y = arange(-philim,philim+step,step)
24   - y = y[:,newaxis]
25   - return exp(-4*log(2) * ((x-center[0])**2 + (y+center[1])**2) / fwhm**2), "jet_"+str(fwhm)
26   -
27   -def computeMap(theta,phi,weight,energy,saveId,nbBins=50,source=isotrop(),borne=[180,180]):
  17 +def computeMap(theta,phi,weight,energy,dir_source,saveId,jet_opening_angle=180,Elim=1e-3,nbBins=50,borne=[180,180]):
28 18 '''
29 19 Plot 2D histogram of arrival direction of the photons with energy upper than Elim
30 20 Input: list of directories
... ... @@ -33,20 +23,22 @@ def computeMap(theta,phi,weight,energy,saveId,nbBins=50,source=isotrop(),borne=[
33 23 Input (optional): limits on theta, phi
34 24 Output: 2D histogram of theta cos(phi), theta sin(phi)
35 25 '''
36   - fig = figure()
37 26  
38   - H,xedges,yedges = histogram2d(theta*cos(phi),theta*sin(phi),bining,weights=weight,range=limits)
39   - conv=convolve2d(H,source[0],boundary="wrap",mode="same")
40   - #H2,xedges,yedges = histogram2d(theta*cos(phi),theta*sin(phi),bining,weights=energy,range=limits)
41   - #conv2=convolve2d(H2,source[0],boundary="wrap",mode="same")
  27 + # apply selection
  28 + cond = (dir_source*degre <= jet_opening_angle) & (energy>Elim) #& (abs(theta)<thetalim) & (abs(phi)<thetalim)
  29 + theta = theta[cond]
  30 + phi = phi[cond]
  31 + weight = weight[cond]
  32 + energy = energy[cond]
  33 +
  34 + fig = figure()
42 35  
43 36 #===============================================================================================#
44 37 # IMAGE
45 38 #===============================================================================================#
46 39 ax = fig.add_subplot(111)
47   - im = ax.matshow(conv,norm=LogNorm(),aspect='auto',
48   - extent=[xedges[0],xedges[-1],yedges[0],yedges[-1]])
49   - #H,xedges,yedges,im = ax.hist2d(theta,phi,bining,weights=weight,norm=LogNorm())
  40 + H,xedges,yedges,im = ax.hist2d(theta*cos(phi),theta*sin(phi),bining,weights=weight,
  41 + norm=LogNorm(),range=limits)
50 42 cbar=fig.colorbar(im, ax=ax)
51 43 cbar.ax.set_ylabel("counts")
52 44 ax.set_xlim((-borne[0],borne[0]))
... ... @@ -54,30 +46,32 @@ def computeMap(theta,phi,weight,energy,saveId,nbBins=50,source=isotrop(),borne=[
54 46 ax.set_xlabel("$\\theta$ $\\cos(\\phi)$ [deg]")
55 47 ax.set_ylabel("$\\theta$ $\\sin(\\phi)$ [deg]")
56 48 ax.grid(b=True,which='major')
57   - ax.set_title(saveId+" - "+source[1])
58   -
59   - savefig(resultDirectory+saveId+"/map_"+source[1]+".png")
  49 + if jet_opening_angle == 180:
  50 + ax.set_title(saveId+" - isotrop")
  51 + savefig(resultDirectory+saveId+"/map_isotrop.png")
  52 + else:
  53 + ax.set_title(saveId+" - jet opening angle: %d degre"%(int(jet_opening_angle)))
  54 + savefig(resultDirectory+saveId+"/map_"+str(int(jet_opening_angle))+".png")
60 55  
  56 + if borne != [180,180]:
  57 + H,xedges,yedges,im = ax.hist2d(theta*cos(phi),theta*sin(phi),bining,weights=weight,
  58 + norm=LogNorm(),range=limits)
61 59 #===============================================================================================#
62 60 # RADIAL DISTRIBUTION AND ENERGY VERSUS ANGLE
63 61 #===============================================================================================#
64   - theta2 = zeros(shape(conv)[0]*shape(conv)[1])
65   - N_evts = zeros(shape(conv)[0]*shape(conv)[1])
66   - #E_evts = zeros(shape(conv)[0]*shape(conv)[1])
  62 + theta2 = zeros(shape(H)[0]*shape(H)[1])
  63 + N_evts = zeros(shape(H)[0]*shape(H)[1])
67 64 theta_cos_phi = xedges[::-1]
68 65 theta_sin_phi = yedges[::-1]
69   - it = nditer(conv,flags=["multi_index"])
70   - #it2 = nditer(conv2,flags=["multi_index"])
  66 + it = nditer(H,flags=["multi_index"])
71 67 k = 0
72 68 while not it.finished:
73 69 i = int(it.multi_index[1])
74 70 j = int(it.multi_index[0])
75 71 theta2[k]= sqrt(theta_cos_phi[i]**2+theta_sin_phi[j]**2)
76 72 N_evts[k]= it[0]
77   - #E_evts[k]= it2[0]
78 73 k += 1
79 74 it.iternext()
80   - #it2.iternext()
81 75  
82 76 theta2=log10(theta2)
83 77 thetamax=max(borne)
... ... @@ -87,18 +81,14 @@ def computeMap(theta,phi,weight,energy,saveId,nbBins=50,source=isotrop(),borne=[
87 81 dtheta = angle[1:nbBins+1]-angle[0:nbBins]
88 82 dndtheta = dn/dtheta*thetacenter
89 83  
90   - #theta2 = 10**theta2
91   - #E_evts = E_evts[E_evts!=0]
92   - #N_evts = N_evts[N_evts!=0]
93   - #theta2 = theta2[E_evts!=0]
94   - #E_evts =log10(E_evts)
95   - #angle,ener =histogram(E_evts,nbBins,weights=N_evts*theta2)
96   - #nb,ener =histogram(E_evts,nbBins,weights=N_evts)
97   - #ener = 10**ener
98   - #enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2
99   - #angle /= nb*degre
  84 + energy =log10(energy)
  85 + angle,ener =histogram(energy,nbBins,weights=weight*theta)
  86 + nb,ener =histogram(energy,nbBins,weights=weight)
  87 + ener = 10**ener
  88 + enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2
  89 + angle /= nb
100 90  
101   - return thetacenter, dndtheta #, enercenter, angle
  91 + return thetacenter, dndtheta , enercenter, angle
102 92  
103 93 def drawRadial(files):
104 94 '''
... ... @@ -121,3 +111,29 @@ def drawRadial(files):
121 111 ax.set_ylabel("$\\theta*dN/d\\theta$")
122 112  
123 113 show()
  114 +
  115 +def drawAngle_vs_energy(files):
  116 + '''
  117 + Plot arrival angle versus energy + analytic expression
  118 + Input: list of directories
  119 + Output: graph of arrival angle versus energy
  120 + '''
  121 + fig = figure()
  122 + ax1 = fig.add_subplot(111)
  123 +
  124 + for fileId in files:
  125 + ener,angle = ReadAngleVsEnergy(fileId,[0,1])
  126 + p=ax1.plot(ener,angle,drawstyle='steps-mid',label=fileId)
  127 +
  128 + yfit = Analytic_theta(ener,fileId)
  129 + ax1.plot(ener,yfit,'--',color="m",linewidth=2,label="analytic")
  130 +
  131 + ax1.legend(loc="best")
  132 + ax1.set_xscale('log')
  133 + ax1.set_yscale('log')
  134 + ax1.set_ylim([0,180])
  135 + ax1.grid(b=True,which='major')
  136 + ax1.set_ylabel("$\\theta$ [rad]")
  137 + ax1.set_xlabel("energy [GeV]")
  138 +
  139 + show()
... ...
Modules/Map.pyc
View file @0b56ded
No preview for this file type
Modules/Spectrum.py
  Diff comments   View file @0b56ded
... ... @@ -32,7 +32,7 @@ def spectrum(energy,weight,generation=[],nbBins=100):
32 32 return enercenter,flux,flux_0
33 33  
34 34  
35   -def drawSpectrum(files,PlotAnalytic=False):
  35 +def drawSpectrum(files,all_source_spectrum=False,PlotAnalytic=False):
36 36 '''
37 37 Plot flux versus energy
38 38 Input: list of directories
... ... @@ -46,25 +46,36 @@ def drawSpectrum(files,PlotAnalytic=False):
46 46 ax2 = fig.add_subplot(gs[1],sharex=ax1)
47 47  
48 48 for fileId in files:
49   - i=1#3
50   - j=1#5
51   - for powerlaw_index in [1,1.5,2,2.5]:
  49 + if all_source_spectrum:
  50 + i=1
  51 + j=1
  52 + for powerlaw_index in [1,1.5,2,2.5]:
  53 + # read files
  54 + Es,Fs = ReadSourceSpectrum(fileId,[0,i])
  55 + p = ax1.plot(Es,Fs,linestyle=':')
  56 + # draw full spectrum
  57 + ener,flux,flux_0 = ReadSpectrum(fileId,[0,j,j+1])
  58 + ax1.plot(ener,flux,color=p[0].get_color(),drawstyle='steps-mid',label="p=%.1f"%powerlaw_index)
  59 + # primary gamma-rays contribution
  60 + ax1.plot(ener,flux_0,color=p[0].get_color(),linestyle='--',drawstyle='steps-mid')
  61 + contrib=flux_0/flux
  62 + ax2.plot(ener,contrib,drawstyle='steps-mid',color=p[0].get_color())
  63 +
  64 + i+=1
  65 + j+=2
  66 +
  67 + else:
52 68 # read files
53   - Es,Fs = ReadSourceSpectrum(fileId,[0,i])
  69 + Es,Fs = ReadSourceSpectrum(fileId,[0,3])
54 70 p = ax1.plot(Es,Fs,linestyle=':')
55 71 # draw full spectrum
56   - ener,flux,flux_0 = ReadSpectrum(fileId,[0,j,j+1])
57   - ax1.plot(ener,flux,color=p[0].get_color(),drawstyle='steps-mid',label="p=%.1f"%powerlaw_index)
  72 + ener,flux,flux_0 = ReadSpectrum(fileId,[0,5,6])
  73 + ax1.plot(ener,flux,color=p[0].get_color(),drawstyle='steps-mid',label=fileId)
58 74 # primary gamma-rays contribution
59 75 ax1.plot(ener,flux_0,color=p[0].get_color(),linestyle='--',drawstyle='steps-mid')
60 76 contrib=flux_0/flux
61 77 ax2.plot(ener,contrib,drawstyle='steps-mid',color=p[0].get_color())
62 78  
63   - i+=1
64   - j+=2
65   -
66   - # add spectrum for Fermi/LAT and CTA (FoV?)
67   -
68 79 if PlotAnalytic:
69 80 minEner=min(ener)
70 81 alpha=flux[ener==minEner]/sqrt(minEner)
... ...
Modules/Spectrum.pyc
View file @0b56ded
No preview for this file type