Thomas Fitoussi / simulation-analysis

Commit b1415bc685f5a8f36338a18fcaaad924d02f2636

Authored by Thomas Fitoussi
1 parent f5a75019
Exists in master

Correction of spectrum and angle distribution 

> all on the same range of energy or angle
Showing 3 changed files with 18 additions and 15 deletions   Show diff stats
Inline Side-by-side
Analysis.py
  Diff comments   View file @b1415bc
... ... @@ -37,6 +37,8 @@ for fileId in argv[1:]:
37 37 phi = phiDir -phiPos
38 38  
39 39 powerlaw_index = 2
  40 + theta_min = 1e-3
  41 + theta_max = 25
40 42  
41 43 #=============================================================================#
42 44 # NO SELECTION
... ... @@ -57,8 +59,7 @@ for fileId in argv[1:]:
57 59 # IMAGING ===================================================================#
58 60 print " ... Computing image"
59 61 nbBins = 50
60   - thetalim = 20
61   - computeMap(theta,phi,weight,energy,fileId,nbBins,borne=[thetalim,thetalim])
  62 + computeMap(theta,phi,weight,energy,fileId,nbBins,borne=[theta_max,theta_max])
62 63  
63 64 # SPECTRUM (MEASURED) =======================================================#
64 65 print " ... Computing spectrum"
... ... @@ -81,7 +82,7 @@ for fileId in argv[1:]:
81 82 # RADIAL ====================================================================#
82 83 print " ... Computing radial distribution"
83 84 nbBins = 100
84   - theta2,dndtheta = radial(theta,weight,nbBins)
  85 + theta2,dndtheta = radial(theta,weight,nbBins,theta_range=[theta_min,theta_max])
85 86 theta2=theta2[:,newaxis]
86 87 dndtheta=dndtheta[:,newaxis]
87 88 Radial = theta2
... ... @@ -90,7 +91,9 @@ for fileId in argv[1:]:
90 91 # TIME DISTRIBUTION AND TIME DELAY VERSUS ANGLE =============================#
91 92 print " ... Computing time distribution"
92 93 nbBins = 100
93   - delta_t,dNdt = timing(time,weight,nbBins)
  94 + min_dt = min(time)
  95 + max_dt = max(time)
  96 + delta_t,dNdt = timing(time,weight,nbBins)#,dt_range=[])
94 97 delta_t=delta_t[:,newaxis]
95 98 dNdt=dNdt[:,newaxis]
96 99 Timing = delta_t
... ... @@ -108,7 +111,7 @@ for fileId in argv[1:]:
108 111  
109 112 # RADIAL =================================================================#
110 113 nbBins = 100
111   - theta2,dndtheta = radial(theta[cond],weight[cond],nbBins)
  114 + theta2,dndtheta = radial(theta[cond],weight[cond],nbBins,theta_range=[theta_min,theta_max])
112 115 dndtheta=dndtheta[:,newaxis]
113 116 Radial = append(Radial,dndtheta,axis=1)
114 117  
... ... @@ -149,7 +152,7 @@ for fileId in argv[1:]:
149 152  
150 153 # RADIAL =================================================================#
151 154 nbBins = 100
152   - theta2,dndtheta = radial(theta[cond],weight[cond],nbBins)
  155 + theta2,dndtheta = radial(theta[cond],weight[cond],nbBins,theta_range=[theta_min,theta_max])
153 156 dndtheta=dndtheta[:,newaxis]
154 157 Radial = append(Radial,dndtheta,axis=1)
155 158  
... ...
Modules/Angle.py
  Diff comments   View file @b1415bc
... ... @@ -5,26 +5,26 @@ from Read import ReadRadial, ReadAngleVsEnergy, ReadGeneration
5 5 from Analytic import Analytic_theta
6 6 from Constants import degre
7 7  
8   -def angle_vs_energy(theta,energy,weight,nbBins=50):
  8 +def angle_vs_energy(theta,energy,weight,nbBins=50,E_range=[1e-3,1e5]):
9 9 '''
10 10 Compute arrival angle versus energy
11 11 Input: directory name
12 12 Output: energy, arrival angle
13 13 '''
14 14 energy =log10(energy)
15   - angle,ener = histogram(energy,nbBins,weights=weight*theta)
16   - nb, ener = histogram(energy,nbBins,weights=weight)
  15 + angle,ener = histogram(energy,nbBins,range=log10(E_range),weights=weight*theta)
  16 + nb, ener = histogram(energy,nbBins,range=log10(E_range),weights=weight)
17 17 ener = 10**ener
18 18 enercenter = (ener[1:nbBins+1]+ener[0:nbBins])/2
19 19 angle /= nb
20 20  
21 21 return enercenter, angle
22 22  
23   -def radial(theta,weight,nbBins=50):
  23 +def radial(theta,weight,nbBins=50,theta_range=[1e-6,180]):
24 24 cond = theta>0
25 25 theta = log10(theta[cond])
26 26 weight = weight[cond]
27   - dn,angle2 = histogram(theta,nbBins,range=[log10(1e-3),log10(25)],weights=weight)
  27 + dn,angle2 = histogram(theta,nbBins,range=log10(theta_range),weights=weight)
28 28 angle2=10**angle2
29 29 #theta = theta**2
30 30 #dn,angle2 = histogram(theta,nbBins,range=[0,25],weights=weight)
... ...
Modules/Spectrum.py
  Diff comments   View file @b1415bc
... ... @@ -5,16 +5,16 @@ from Read import ReadSpectrum, ReadSourceSpectrum, ReadGeneration
5 5 from Constants import degre
6 6 from Analytic import Ethreshold_gg, Analytic_flux, Eic
7 7  
8   -def spectrum(energy,weight,nbBins=100):
  8 +def spectrum(energy,weight,nbBins=100,E_range=[1e-3,1e5]):
9 9 '''
10 10 Compute flux versus energy
11 11 Input: events energy and weight (optional: number of bins)
12 12 Output: energy, flux (E^2 dN/dE)
13 13 '''
14   - Emin=min(energy)
15   - Emax=max(energy)
  14 + Emin=E_range[0]#min(energy)
  15 + Emax=E_range[1]#max(energy)
16 16 energy =log10(energy)
17   - flux,ener=histogram(energy,nbBins,range=[log10(Emin),log10(Emax)],weights=weight)
  17 + flux,ener=histogram(energy,nbBins,range=log10(E_range),weights=weight)
18 18 ener = 10**ener
19 19 enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2
20 20 binSize=ener[1:nbBins+1]-ener[0:nbBins]
... ...