diff --git a/Analysis.py b/Analysis.py index 59167a1..e47b358 100755 --- a/Analysis.py +++ b/Analysis.py @@ -37,6 +37,8 @@ for fileId in argv[1:]: phi = phiDir -phiPos powerlaw_index = 2 + theta_min = 1e-3 + theta_max = 25 #=============================================================================# # NO SELECTION @@ -57,8 +59,7 @@ for fileId in argv[1:]: # IMAGING ===================================================================# print " ... Computing image" nbBins = 50 - thetalim = 20 - computeMap(theta,phi,weight,energy,fileId,nbBins,borne=[thetalim,thetalim]) + computeMap(theta,phi,weight,energy,fileId,nbBins,borne=[theta_max,theta_max]) # SPECTRUM (MEASURED) =======================================================# print " ... Computing spectrum" @@ -81,7 +82,7 @@ for fileId in argv[1:]: # RADIAL ====================================================================# print " ... Computing radial distribution" nbBins = 100 - theta2,dndtheta = radial(theta,weight,nbBins) + theta2,dndtheta = radial(theta,weight,nbBins,theta_range=[theta_min,theta_max]) theta2=theta2[:,newaxis] dndtheta=dndtheta[:,newaxis] Radial = theta2 @@ -90,7 +91,9 @@ for fileId in argv[1:]: # TIME DISTRIBUTION AND TIME DELAY VERSUS ANGLE =============================# print " ... Computing time distribution" nbBins = 100 - delta_t,dNdt = timing(time,weight,nbBins) + min_dt = min(time) + max_dt = max(time) + delta_t,dNdt = timing(time,weight,nbBins)#,dt_range=[]) delta_t=delta_t[:,newaxis] dNdt=dNdt[:,newaxis] Timing = delta_t @@ -108,7 +111,7 @@ for fileId in argv[1:]: # RADIAL =================================================================# nbBins = 100 - theta2,dndtheta = radial(theta[cond],weight[cond],nbBins) + theta2,dndtheta = radial(theta[cond],weight[cond],nbBins,theta_range=[theta_min,theta_max]) dndtheta=dndtheta[:,newaxis] Radial = append(Radial,dndtheta,axis=1) @@ -149,7 +152,7 @@ for fileId in argv[1:]: # RADIAL =================================================================# nbBins = 100 - theta2,dndtheta = radial(theta[cond],weight[cond],nbBins) + theta2,dndtheta = radial(theta[cond],weight[cond],nbBins,theta_range=[theta_min,theta_max]) dndtheta=dndtheta[:,newaxis] Radial = append(Radial,dndtheta,axis=1) diff --git a/Modules/Angle.py b/Modules/Angle.py index c38d7f6..b4d70b8 100644 --- a/Modules/Angle.py +++ b/Modules/Angle.py @@ -5,26 +5,26 @@ from Read import ReadRadial, ReadAngleVsEnergy, ReadGeneration from Analytic import Analytic_theta from Constants import degre -def angle_vs_energy(theta,energy,weight,nbBins=50): +def angle_vs_energy(theta,energy,weight,nbBins=50,E_range=[1e-3,1e5]): ''' Compute arrival angle versus energy Input: directory name Output: energy, arrival angle ''' energy =log10(energy) - angle,ener = histogram(energy,nbBins,weights=weight*theta) - nb, ener = histogram(energy,nbBins,weights=weight) + angle,ener = histogram(energy,nbBins,range=log10(E_range),weights=weight*theta) + nb, ener = histogram(energy,nbBins,range=log10(E_range),weights=weight) ener = 10**ener enercenter = (ener[1:nbBins+1]+ener[0:nbBins])/2 angle /= nb return enercenter, angle -def radial(theta,weight,nbBins=50): +def radial(theta,weight,nbBins=50,theta_range=[1e-6,180]): cond = theta>0 theta = log10(theta[cond]) weight = weight[cond] - dn,angle2 = histogram(theta,nbBins,range=[log10(1e-3),log10(25)],weights=weight) + dn,angle2 = histogram(theta,nbBins,range=log10(theta_range),weights=weight) angle2=10**angle2 #theta = theta**2 #dn,angle2 = histogram(theta,nbBins,range=[0,25],weights=weight) diff --git a/Modules/Spectrum.py b/Modules/Spectrum.py index cb651ec..2f785e0 100644 --- a/Modules/Spectrum.py +++ b/Modules/Spectrum.py @@ -5,16 +5,16 @@ from Read import ReadSpectrum, ReadSourceSpectrum, ReadGeneration from Constants import degre from Analytic import Ethreshold_gg, Analytic_flux, Eic -def spectrum(energy,weight,nbBins=100): +def spectrum(energy,weight,nbBins=100,E_range=[1e-3,1e5]): ''' Compute flux versus energy Input: events energy and weight (optional: number of bins) Output: energy, flux (E^2 dN/dE) ''' - Emin=min(energy) - Emax=max(energy) + Emin=E_range[0]#min(energy) + Emax=E_range[1]#max(energy) energy =log10(energy) - flux,ener=histogram(energy,nbBins,range=[log10(Emin),log10(Emax)],weights=weight) + flux,ener=histogram(energy,nbBins,range=log10(E_range),weights=weight) ener = 10**ener enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2 binSize=ener[1:nbBins+1]-ener[0:nbBins] -- libgit2 0.21.2