From f5a7501927b775d30ec58a0d27b8b508657e8aec Mon Sep 17 00:00:00 2001 From: Thomas Fitoussi Date: Mon, 12 Oct 2015 17:19:35 +0200 Subject: [PATCH] Rewrite Analysis.py : > only one power law spectrum index > compute results gen by gen and by energy range when it makes sense > add observables.py to plot 3D space of arrival parameters (E, dt, theta) --- .gitignore | 2 ++ Analysis.py | 296 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- EBL-spectrums.py | 4 ++-- Modules/Analytic.pyc | Bin 5449 -> 0 bytes Modules/Angle.py | 56 ++++++++++++++++++++++++++++---------------------------- Modules/Angle.pyc | Bin 3897 -> 0 bytes Modules/Constants.pyc | Bin 1198 -> 0 bytes Modules/Generation.pyc | Bin 1175 -> 0 bytes Modules/Integrand.pyc | Bin 1067 -> 0 bytes Modules/Map.py | 10 +++------- Modules/Map.pyc | Bin 2704 -> 0 bytes Modules/Observables.py | 56 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Modules/Read.py | 2 +- Modules/Read.pyc | Bin 3938 -> 0 bytes Modules/Spectrum.py | 96 +++++++++++++++++++++++++++++++----------------------------------------------------------------- Modules/Spectrum.pyc | Bin 3420 -> 0 bytes Modules/Timing.py | 91 ++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------- Modules/Timing.pyc | Bin 5673 -> 0 bytes Modules/Weight.pyc | Bin 2434 -> 0 bytes Modules/__init__.py | 1 + Modules/__init__.pyc | Bin 110 -> 0 bytes Results/.gitignore | 2 ++ lambda_gg.py | 57 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 23 files changed, 327 insertions(+), 346 deletions(-) create mode 100644 .gitignore delete mode 100644 Modules/Analytic.pyc delete mode 100644 Modules/Angle.pyc delete mode 100644 Modules/Constants.pyc delete mode 100644 Modules/Generation.pyc delete mode 100644 Modules/Integrand.pyc delete mode 100644 Modules/Map.pyc create mode 100644 Modules/Observables.py delete mode 100644 Modules/Read.pyc delete mode 100644 Modules/Spectrum.pyc delete mode 100644 Modules/Timing.pyc delete mode 100644 Modules/Weight.pyc delete mode 100644 Modules/__init__.pyc create mode 100644 Results/.gitignore create mode 100755 lambda_gg.py diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..c9b568f --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +*.pyc +*.swp diff --git a/Analysis.py b/Analysis.py index b705273..59167a1 100755 --- a/Analysis.py +++ b/Analysis.py @@ -10,21 +10,18 @@ from Modules.Angle import angle_vs_energy, radial from Modules.Timing import timing from Modules.Constants import degre -def InProgress(i,Nmax): - print " ", (i*100)/Nmax, "% done" - if shape(argv)[0] < 2: print "not enough arguments (at least 1)" exit() -PowerSpectrum=[1,1.5,2,2.5] +powerlaw_index=2 for fileId in argv[1:]: - print "#=============================================================================#" - print "# Analysis of", fileId - print "#=============================================================================#" + print "#===============================================#" + print "# Analysis of", fileId + print "#===============================================#" # read files - print "# 1. Reading data" + print " > Reading data" time = ReadTime(fileId) energy = ReadEnergy(fileId) weightini, generation, theta_arrival, Esource, dir_source = ReadExtraFile(fileId,[2,3,4,5,6]) @@ -39,228 +36,135 @@ for fileId in argv[1:]: theta = thetaDir - thetaPos phi = phiDir -phiPos - Gen_contrib = [] - Gen_cont = zeros((int(max(generation))+1)) - Source = [] - Spectrum = [] - Timing = [] - Angle_Energy = [] - Radial = [] - - print "# 2. Computing powerlaw spectrum" - #=============================================================================# - # VERSUS SOURCE SPECTRUM - #=============================================================================# - for powerlaw_index in PowerSpectrum: - # apply source spectrum - weight_source = (Esource/min(Esource))**(1-powerlaw_index) - weight = weightini* weight_source - - # GENERATION CONTRIBUTION ====================================================# - NbTotEvents = sum(weight) - if Gen_contrib==[]: - Gen_contrib = arange(0,max(generation)+1,1)[:,newaxis] - - for gen in list(set(generation)): - Gen_cont[int(gen)] = sum(weight[generation==gen])/NbTotEvents *100 - Gen_contrib = append(Gen_contrib,Gen_cont[:,newaxis],axis=1) - - # SPECTRUM (SOURCE AND MEASURED) ============================================# - nbBins = 100 - # draw source spectrum - Es=array(list(set(Esource))) - Ws= (Es/min(Es))**(1-powerlaw_index) / ratio - Es,Fs = spectrum(Es,Ws,nbBins=nbBins) - Es=Es[:,newaxis] - Fs=Fs[:,newaxis] - if Source==[]: - Source = Es - Source = append(Source,Fs,axis=1) - else: - Source = append(Source,Fs,axis=1) - - # primary gamma-rays contribution and full spectrum - ener,flux,flux_0 = spectrum(energy,weight,generation,nbBins) - ener=ener[:,newaxis] - flux=flux[:,newaxis] - flux_0=flux_0[:,newaxis] - if Spectrum==[]: - Spectrum = ener - Spectrum = append(Spectrum,flux,axis=1) - Spectrum = append(Spectrum,flux_0,axis=1) - else: - Spectrum = append(Spectrum,flux,axis=1) - Spectrum = append(Spectrum,flux_0,axis=1) - - # IMAGING, RADIAL DISTRIBUTION AND ANGLE VERSUS ENERGY ======================# - nbBins = 50 - cond = energy>1e-3 - ener,angle = angle_vs_energy(theta_arrival[cond],energy[cond],weight[cond],nbBins) - theta2,dndtheta2 = radial(theta[cond],weight[cond],nbBins) - theta2=theta2[:,newaxis] - dndtheta2=dndtheta2[:,newaxis] - ener=ener[:,newaxis] - angle=angle[:,newaxis] - - if Radial==[]: - Radial = theta2 - Radial = append(Radial,dndtheta2,axis=1) - Angle_Energy = ener - Angle_Energy = append(Angle_Energy,angle,axis=1) - else: - Radial = append(Radial,dndtheta2,axis=1) - Angle_Energy = append(Angle_Energy,angle,axis=1) - - - # TIME DISTRIBUTION AND TIME DELAY VERSUS ANGLE =============================# - nbBins = 100 - delta_t,dNdt = timing(time,weight,nbBins) - delta_t=delta_t[:,newaxis] - dNdt=dNdt[:,newaxis] - if Timing==[]: - Timing = delta_t - Timing = append(Timing,dNdt,axis=1) - else: - Timing = append(Timing,dNdt,axis=1) - - InProgress(PowerSpectrum.index(powerlaw_index)+1,shape(PowerSpectrum)[0]) + powerlaw_index = 2 #=============================================================================# - # VERSUS JET OPENING ANGLE + # NO SELECTION #=============================================================================# - Jet_Opening=[180,60,30] # degre (180 <=> isotrop) - powerlaw_index = 2 - Elim = 1e-1 # GeV - thetalim = 20 # degre - - # apply source spectrum + print " > No selection of events ..." + # SPECTRUM (SOURCE) =========================================================# weight_source = (Esource/min(Esource))**(1-powerlaw_index) weight = weightini* weight_source - - print "# 3. compute images and radial distribution" - for jet_opening_angle in Jet_Opening: - # apply selection ( /!\ USE DECREASING VALUES OF jet_opening_angle ) - cond = (dir_source*degre <= jet_opening_angle) #& (energy>Elim) #& (abs(theta) Selection by energy range ..." Emin = [1e-3,1e0,1e3] #GeV Emax = [1e0,1e3,1e5] #GeV + for n in arange(0,3,1): cond= (energy>Emin[n]) & (energy Selection by generation ..." + NbTotEvents = sum(weight) gen_tab =list(set(generation)) - powerlaw_index = 2 - weight_source = (Esource/min(Esource))**(1-powerlaw_index) - weight = weightini* weight_source - + Gen_cont = zeros([shape(gen_tab)[0],2]) for gen in gen_tab: cond = generation==gen - print shape(time[cond]) + contrib =sum(weight[cond])/NbTotEvents*100 + i = gen_tab.index(gen) + Gen_cont[i,0]=int(gen) + Gen_cont[i,1]=contrib + print " ... gen=",int(gen),"-> contribution:",int(contrib),"%" + + # SPECTRUM (MEASURED) ====================================================# + nbBins = 100 + ener,flux = spectrum(energy[cond],weight[cond],nbBins) + flux=flux[:,newaxis] + Spectrum = append(Spectrum,flux,axis=1) - # TIME DISTRIBUTION AND TIME DELAY VERSUS ANGLE =============================# + # ANGLE VERSUS ENERGY ====================================================# + nbBins = 100 + ener,angle = angle_vs_energy(theta_arrival[cond],energy[cond],weight[cond],nbBins) + angle=angle[:,newaxis] + Angle_Energy = append(Angle_Energy,angle,axis=1) + + # RADIAL =================================================================# + nbBins = 100 + theta2,dndtheta = radial(theta[cond],weight[cond],nbBins) + dndtheta=dndtheta[:,newaxis] + Radial = append(Radial,dndtheta,axis=1) + + # TIME DISTRIBUTION AND TIME DELAY VERSUS ANGLE ==========================# nbBins = 100 delta_t,dNdt = timing(time[cond],weight[cond],nbBins) - delta_t=delta_t[:,newaxis] dNdt=dNdt[:,newaxis] Timing = append(Timing,dNdt,axis=1) #=============================================================================# - print "# 4. writing files" - savetxt(resultDirectory+fileId+"/Generation.txt",Gen_contrib) + print " > writing files" savetxt(resultDirectory+fileId+"/Spectrum.txt",Spectrum) savetxt(resultDirectory+fileId+"/Source_spectrum.txt",Source) savetxt(resultDirectory+fileId+"/Angle_vs_Energy.txt",Angle_Energy) savetxt(resultDirectory+fileId+"/Radial_distribution.txt",Radial) savetxt(resultDirectory+fileId+"/Timing.txt",Timing) - print "#=============================================================================#" - + savetxt(resultDirectory+fileId+"/Generation.txt",Gen_cont) + print "#===============================================#" diff --git a/EBL-spectrums.py b/EBL-spectrums.py index 603bdfc..3e36123 100755 --- a/EBL-spectrums.py +++ b/EBL-spectrums.py @@ -58,7 +58,7 @@ ax.plot(hv,nCMB(hv,z)*hv**2,"-k",label="CMB") ax.set_xscale('log') ax.set_yscale('log') ax.grid(b=True,which='major') -ax.legend(loc="best",frameon=False,framealpha=0.5) +ax.legend(loc="best")#,frameon=False,framealpha=0.5) ax.set_xlabel("energy [eV]") ax.set_ylabel("$n$ [photon.eV.cm$^{-3}$]") @@ -80,7 +80,7 @@ for z_index in [1,9,12,15,17]: ax2.set_xscale('log') ax2.set_yscale('log') ax2.grid(b=True,which='major') -ax2.legend(loc="best",frameon=False,framealpha=0.5) +ax2.legend(loc="best")#,frameon=False,framealpha=0.5) ax2.set_xlabel("energy [eV]") ax2.set_ylabel("$n$ [photon.eV.cm$^{-3}$]") diff --git a/Modules/Analytic.pyc b/Modules/Analytic.pyc deleted file mode 100644 index 2bd40ab..0000000 Binary files a/Modules/Analytic.pyc and /dev/null differ diff --git a/Modules/Angle.py b/Modules/Angle.py index 5d52a49..c38d7f6 100644 --- a/Modules/Angle.py +++ b/Modules/Angle.py @@ -1,7 +1,7 @@ -from numpy import log10, histogram +from numpy import log10, histogram, arange from matplotlib.pyplot import figure, show from matplotlib import gridspec -from Read import ReadRadial, ReadAngleVsEnergy +from Read import ReadRadial, ReadAngleVsEnergy, ReadGeneration from Analytic import Analytic_theta from Constants import degre @@ -21,7 +21,7 @@ def angle_vs_energy(theta,energy,weight,nbBins=50): return enercenter, angle def radial(theta,weight,nbBins=50): - cond = theta!=0 + cond = theta>0 theta = log10(theta[cond]) weight = weight[cond] dn,angle2 = histogram(theta,nbBins,range=[log10(1e-3),log10(25)],weights=weight) @@ -35,7 +35,7 @@ def radial(theta,weight,nbBins=50): return thetacenter, dndtheta -def drawRadial(files,all_source_spectrum=False,all_jets=False): +def drawRadial(files,plot_gen=False,plot_energy_range=False): ''' Plot flux versus arrival angle Input: list of directories @@ -45,21 +45,24 @@ def drawRadial(files,all_source_spectrum=False,all_jets=False): ax = fig.add_subplot(111) for fileId in files: - if all_source_spectrum: - i=1 - for powerlaw_index in [1,1.5,2,2.5]: + theta,dndtheta2 = ReadRadial(fileId,[0,1]) + ax.plot(theta,dndtheta2,drawstyle='steps-mid') + + if plot_energy_range: + labels = ["MeV band","GeV band","TeV band"] + i = 2 + for n in arange(0,3,1): theta,dndtheta2 = ReadRadial(fileId,[0,i]) - ax.plot(theta,dndtheta2,drawstyle='steps-mid',label="p=%.1f"%powerlaw_index) + ax.plot(theta,dndtheta2,drawstyle='steps-mid',label=labels[n]) i+=1 - elif all_jets: + + if plot_gen: + generation = ReadGeneration(fileId,[0]) i=5 - for jet_opening in ["isotrop","60 deg","30 deg"]: + for gen in generation: theta,dndtheta2 = ReadRadial(fileId,[0,i]) - ax.plot(theta,dndtheta2,drawstyle='steps-mid',label=jet_opening) + ax.plot(theta,dndtheta2,drawstyle='steps-mid',label="gen=%.0f"%gen) i+=1 - else: - theta,dndtheta2 = ReadRadial(fileId,[0,5]) - ax.plot(theta,dndtheta2,drawstyle='steps-mid',label=fileId) ax.legend(loc="best") ax.set_xscale('log') @@ -70,7 +73,7 @@ def drawRadial(files,all_source_spectrum=False,all_jets=False): show() -def drawAngle_vs_energy(files,all_source_spectrum=False,all_jets=False,PlotAnalytic=False): +def drawAngle_vs_energy(files,plot_gen=False,plot_analytic=False): ''' Plot arrival angle versus energy + analytic expression Input: list of directories @@ -80,21 +83,18 @@ def drawAngle_vs_energy(files,all_source_spectrum=False,all_jets=False,PlotAnaly ax1 = fig.add_subplot(111) for fileId in files: - if all_source_spectrum: - i=1 - for powerlaw_index in [1,1.5,2,2.5]: - ener,angle = ReadAngleVsEnergy(fileId,[0,i]) - p=ax1.plot(ener,angle,drawstyle='steps-mid',label="p=%.1f"%powerlaw_index) - elif all_jets: - i=5 - for jet_opening in ["isotrop","60 deg","30 deg"]: + ener,angle = ReadAngleVsEnergy(fileId,[0,1]) + p=ax1.plot(ener,angle,drawstyle='steps-mid') + + if plot_gen: + generation = ReadGeneration(fileId,[0]) + i=2 + for gen in generation: ener,angle = ReadAngleVsEnergy(fileId,[0,i]) - p=ax1.plot(ener,angle,drawstyle='steps-mid',label=jet_opening) - else: - ener,angle = ReadAngleVsEnergy(fileId,[0,1]) - p=ax1.plot(ener,angle,drawstyle='steps-mid',label=fileId) + ax1.plot(ener,angle,drawstyle='steps-mid',label="gen=%.0f"%gen) + i+=1 - if PlotAnalytic: + if plot_analytic: yfit = Analytic_theta(ener,fileId) ax1.plot(ener,yfit,'--',color=p[0].get_color(),linewidth=2) diff --git a/Modules/Angle.pyc b/Modules/Angle.pyc deleted file mode 100644 index 9edcd34..0000000 Binary files a/Modules/Angle.pyc and /dev/null differ diff --git a/Modules/Constants.pyc b/Modules/Constants.pyc deleted file mode 100644 index 3bd6752..0000000 Binary files a/Modules/Constants.pyc and /dev/null differ diff --git a/Modules/Generation.pyc b/Modules/Generation.pyc deleted file mode 100644 index b71b40b..0000000 Binary files a/Modules/Generation.pyc and /dev/null differ diff --git a/Modules/Integrand.pyc b/Modules/Integrand.pyc deleted file mode 100644 index b31c093..0000000 Binary files a/Modules/Integrand.pyc and /dev/null differ diff --git a/Modules/Map.py b/Modules/Map.py index 9c7cf6b..f0b46e7 100644 --- a/Modules/Map.py +++ b/Modules/Map.py @@ -15,7 +15,7 @@ philim=float(philim) thetalim=float(thetalim) limits = [[-thetalim,thetalim],[-philim,philim]] -def computeMap(theta,phi,weight,energy,saveId,jet_opening_angle=180,nbBins=50,borne=[180,180]): +def computeMap(theta,phi,weight,energy,saveId,nbBins=50,borne=[180,180]): ''' Plot 2D histogram of arrival direction of the photons with energy upper than Elim Input: list of directories @@ -41,12 +41,8 @@ def computeMap(theta,phi,weight,energy,saveId,jet_opening_angle=180,nbBins=50,bo ax.set_xlabel("$\\theta$ $\\cos(\\phi)$ [deg]") ax.set_ylabel("$\\theta$ $\\sin(\\phi)$ [deg]") ax.grid(b=True,which='major') - if jet_opening_angle == 180: - ax.set_title(saveId+" - isotrop") - savefig(resultDirectory+saveId+"/map_isotrop.png") - else: - ax.set_title(saveId+" - jet opening angle: %d degre"%(int(jet_opening_angle))) - savefig(resultDirectory+saveId+"/map_"+str(int(jet_opening_angle))+".png") + ax.set_title(saveId+" - isotrop") + savefig(resultDirectory+saveId+"/map_isotrop.png") #def isotrop(): # return array([[1]]), "isotropic" diff --git a/Modules/Map.pyc b/Modules/Map.pyc deleted file mode 100644 index 7388a80..0000000 Binary files a/Modules/Map.pyc and /dev/null differ diff --git a/Modules/Observables.py b/Modules/Observables.py new file mode 100644 index 0000000..5b26d8d --- /dev/null +++ b/Modules/Observables.py @@ -0,0 +1,56 @@ +from numpy import sort, log10 +from matplotlib.pyplot import figure, show +from mpl_toolkits.mplot3d import Axes3D +from Modules.Read import ReadEnergy, ReadTime, ReadExtraFile +from Modules.Constants import degre + +def drawObservablesSpace(fileId): + fig = figure() + ax = fig.add_subplot(232,projection='3d') + ax1 = fig.add_subplot(233) + ax2 = fig.add_subplot(231) + ax3 = fig.add_subplot(235) + + energy = ReadEnergy(fileId) + delay = ReadTime(fileId) + generation, theta = ReadExtraFile(fileId,[3,4]) + theta*=degre + + colors=['b','g','r','m','c','y'] + i=0 + + for gen in sort(list(set(generation))): + cond = (generation==gen) & (delay>0) & (theta>0) + ax.scatter(log10(energy[cond]),log10(theta[cond]),log10(delay[cond]),color=colors[i],label="gen = %.0f"%gen) + ax1.scatter(energy[cond],delay[cond],color=colors[i]) + ax2.scatter(theta[cond],delay[cond],color=colors[i]) + ax3.scatter(theta[cond],energy[cond],color=colors[i]) + i+=1 + + ax.legend(loc="best") + #ax.set_xscale('log') + #ax.set_yscale('log') + #ax.set_zscale('log') + ax.set_xlabel("energy [GeV]") + ax.set_ylabel("$\\theta$ [deg]") + ax.set_zlabel("Time delay [s]") + + ax1.grid(b=True,which='major') + ax1.set_xscale('log') + ax1.set_yscale('log') + ax1.set_xlabel("energy [GeV]") + ax1.set_ylabel("Time delay [s]") + + ax2.grid(b=True,which='major') + ax2.set_xscale('log') + ax2.set_yscale('log') + ax2.set_xlabel("$\\theta$ [deg]") + ax2.set_ylabel("Time delay [s]") + + ax3.grid(b=True,which='major') + ax3.set_xscale('log') + ax3.set_yscale('log') + ax3.set_xlabel("$\\theta$ [deg]") + ax3.set_ylabel("energy [GeV]") + + show() diff --git a/Modules/Read.py b/Modules/Read.py index a43c725..d89ee76 100644 --- a/Modules/Read.py +++ b/Modules/Read.py @@ -60,5 +60,5 @@ def ReadRadial(fileName,cols=[0,1,2,3,4]): def ReadTiming(fileName,cols=[0,1,2,3,4]): return loadtxt(resultDirectory+fileName+"/Timing.txt",unpack=True,usecols=cols) -def ReadGeneration(fileName,cols=[0,1,2,3,4]): +def ReadGeneration(fileName,cols=[0,1]): return loadtxt(resultDirectory+fileName+"/Generation.txt",unpack=True,usecols=cols) diff --git a/Modules/Read.pyc b/Modules/Read.pyc deleted file mode 100644 index e4d6b42..0000000 Binary files a/Modules/Read.pyc and /dev/null differ diff --git a/Modules/Spectrum.py b/Modules/Spectrum.py index ef42fd0..cb651ec 100644 --- a/Modules/Spectrum.py +++ b/Modules/Spectrum.py @@ -1,11 +1,11 @@ from numpy import histogram, log10, sqrt, array from matplotlib.pyplot import figure, show, setp from matplotlib import gridspec -from Read import ReadSpectrum, ReadSourceSpectrum +from Read import ReadSpectrum, ReadSourceSpectrum, ReadGeneration from Constants import degre -from Analytic import Ethreshold_gg +from Analytic import Ethreshold_gg, Analytic_flux, Eic -def spectrum(energy,weight,generation=[],nbBins=100): +def spectrum(energy,weight,nbBins=100): ''' Compute flux versus energy Input: events energy and weight (optional: number of bins) @@ -19,20 +19,9 @@ def spectrum(energy,weight,generation=[],nbBins=100): enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2 binSize=ener[1:nbBins+1]-ener[0:nbBins] flux = enercenter**2 * flux/binSize + return enercenter,flux - if generation == []: - return enercenter,flux - else: - cond=(generation==0) - flux_0,ener_0=histogram(energy[cond],nbBins,range=[log10(Emin),log10(Emax)],weights=weight[cond]) - ener_0 = 10**ener_0 - enercenter_0=(ener_0[1:nbBins+1]+ener_0[0:nbBins])/2 - binSize=ener_0[1:nbBins+1]-ener_0[0:nbBins] - flux_0 = enercenter_0**2 * flux_0/binSize - return enercenter,flux,flux_0 - - -def drawSpectrum(files,all_source_spectrum=False,all_jets=False,PlotAnalytic=False): +def drawSpectrum(files,plot_gen=False,plot_source=False,plot_analytic=False): ''' Plot flux versus energy Input: list of directories @@ -47,59 +36,36 @@ def drawSpectrum(files,all_source_spectrum=False,all_jets=False,PlotAnalytic=Fal ax1 = fig.add_subplot(111) for fileId in files: - if all_source_spectrum: - i=1 - j=1 - for powerlaw_index in [1,1.5,2,2.5]: - # read files - Es,Fs = ReadSourceSpectrum(fileId,[0,i]) - p = ax1.plot(Es,Fs,linestyle=':') - # draw full spectrum - ener,flux,flux_0 = ReadSpectrum(fileId,[0,j,j+1]) - ax1.plot(ener,flux,color=p[0].get_color(),drawstyle='steps-mid',label="p=%.1f"%powerlaw_index) - # primary gamma-rays contribution - #ax1.plot(ener,flux_0,color=p[0].get_color(),linestyle='--',drawstyle='steps-mid') - #contrib=flux_0/flux - #ax2.plot(ener,contrib,drawstyle='steps-mid',color=p[0].get_color()) - i+=1 - j+=2 - elif all_jets: - i=5 - j=9 - for jet_opening in ["isotrop","60 deg","30 deg"]: + # draw full spectrum + ener,flux = ReadSpectrum(fileId,[0,1]) + p = ax1.plot(ener,flux,drawstyle='steps-mid') + + if plot_gen: + generation = ReadGeneration(fileId,[0]) + j=2 + for gen in generation: # read files - Es,Fs = ReadSourceSpectrum(fileId,[0,i]) - p = ax1.plot(Es,Fs,linestyle=':') - # draw full spectrum - ener,flux,flux_0 = ReadSpectrum(fileId,[0,j,j+1]) - ax1.plot(ener,flux,color=p[0].get_color(),drawstyle='steps-mid',label=jet_opening) - # primary gamma-rays contribution - #ax1.plot(ener,flux_0,color=p[0].get_color(),linestyle='--',drawstyle='steps-mid') - #contrib=flux_0/flux - #ax2.plot(ener,contrib,drawstyle='steps-mid',color=p[0].get_color()) - i+=1 - j+=2 + ener,flux = ReadSpectrum(fileId,[0,j]) + ax1.plot(ener,flux,drawstyle='steps-mid',label="gen=%.0f"%gen) + j+=1 - else: - # read files - #Es,Fs = ReadSourceSpectrum(fileId,[0,3]) - #p = ax1.plot(Es,Fs,linestyle=':') - # draw full spectrum - ener,flux,flux_0 = ReadSpectrum(fileId,[0,5,6]) - ax1.plot(ener,flux,drawstyle='steps-mid') - #ax1.plot(ener,flux,color=p[0].get_color(),drawstyle='steps-mid',label=fileId) - # primary gamma-rays contribution - #ax1.plot(ener,flux_0,color=p[0].get_color(),linestyle='--',drawstyle='steps-mid') - #contrib=flux_0/flux - #ax2.plot(ener,contrib,drawstyle='steps-mid',color=p[0].get_color()) + if plot_source: + Es,Fs = ReadSourceSpectrum(fileId,[0,1]) + ax1.plot(Es,Fs,color=p[0].get_color(),linestyle=':') - if PlotAnalytic: + if plot_analytic: minEner=min(ener) - alpha=flux[ener==minEner]/sqrt(minEner) - ax1.plot(ener,alpha*sqrt(ener),'--m',linewidth=2) - ax1.plot(ener,ener**(0)*max(flux)*0.95,'--m',linewidth=2) - ax1.axvline(x=Ethreshold_gg(), ymin=0., ymax = 1., color='r', linewidth=2) - ax1.axvline(x=1e4, ymin=0., ymax = 1., color='g', linewidth=2) + #alpha=flux[ener==minEner]/sqrt(minEner) + #ax1.plot(ener,alpha*sqrt(ener/100),'.-g',linewidth=2) + ax1.plot(ener,2*Analytic_flux(ener),'--g',linewidth=1) + Elim= Eic(Ethreshold_gg()/2) + alpha = 2*4.63e3/Elim**(1/4.) + ax1.plot(ener,alpha*sqrt(ener),'--r',linewidth=1) + alpha = 2*14.6e3 + ax1.plot(ener,alpha*(ener/100)**(1/4.),'--m',linewidth=1) + #ax1.plot(ener,ener**(0)*max(flux)*0.95,'--m',linewidth=2) + #ax1.axvline(x=Ethreshold_gg(), ymin=0., ymax = 1., color='r', linewidth=2) + #ax1.axvline(x=1e4, ymin=0., ymax = 1., color='g', linewidth=2) ax1.legend(loc="best") ax1.set_xscale('log') diff --git a/Modules/Spectrum.pyc b/Modules/Spectrum.pyc deleted file mode 100644 index 2593ae4..0000000 Binary files a/Modules/Spectrum.pyc and /dev/null differ diff --git a/Modules/Timing.py b/Modules/Timing.py index d34759b..852e4ba 100644 --- a/Modules/Timing.py +++ b/Modules/Timing.py @@ -1,6 +1,6 @@ -from numpy import histogram, log10, shape, arange, pi, convolve +from numpy import histogram, log10, shape, arange, pi, convolve, r_ from matplotlib.pyplot import figure, show -from Read import ReadTime, ReadTiming, ReadExtraFile, ReadProfile, ReadEnergy +from Read import ReadTime, ReadTiming, ReadExtraFile, ReadProfile, ReadEnergy, ReadGeneration from Constants import yr, degre from Analytic import Analytic_delay_vs_theta @@ -23,50 +23,35 @@ def timing(time,weight,nbBins=100): binSize=dt[1:nbBins+1]-dt[0:nbBins] dNdt=dN/binSize *timecenter - return timecenter, dNdt + return timecenter, dNdt#/max(dNdt) -def drawTiming(files,plot_case=0): +def drawTiming(files,plot_gen=False,plot_energy_range=False): ''' Plot flux versus time delay Input: list of directories - plot_case: (0) nothing, (1) all source spectrum, (2) all jets, (3) by energy band, - (4) by generation Output: graph of flux versus time delay ''' fig = figure() ax = fig.add_subplot(111) for fileId in files: - if plot_case == 1: - i=1 - for powerlaw_index in [1,1.5,2,2.5]: - delta_t,dNdt = ReadTiming(fileId,[0,i]) - ax.plot(delta_t,dNdt,linestyle="steps-mid",label="p=%.1f"%powerlaw_index) - i+=1 - elif plot_case == 2: - i=5 - for jet_opening in ["isotrop","60 degre","30 degre"]: - delta_t,dNdt = ReadTiming(fileId,[0,i]) - ax.plot(delta_t,dNdt,linestyle="steps-mid",label=jet_opening) - i+=1 - elif plot_case == 3: - i=8 - Emin = [1e-3,1e0,1e3] #GeV - Emax = [1e0,1e3,1e5] #GeV + delta_t,dNdt = ReadTiming(fileId,[0,1]) + p=ax.plot(delta_t,dNdt,linestyle="steps-mid") + + if plot_energy_range: + labels = ["MeV band","GeV band","TeV band"] + i = 2 for n in arange(0,3,1): - label = "%1.0e Gev"%Emin[n]+"< $E_\\gamma$ < %1.0e Gev"%Emax[n] delta_t,dNdt = ReadTiming(fileId,[0,i]) - ax.plot(delta_t,dNdt,linestyle="steps-mid",label=label) + ax.plot(delta_t,dNdt,linestyle="steps-mid",label=labels[n]) i+=1 - elif plot_case == 4: - i=11 - generation = ReadExtraFile(fileId,[3]) - for gen in list(set(generation)): + + if plot_gen: + generation = ReadGeneration(fileId,[0]) + i=5 + for gen in generation: delta_t,dNdt = ReadTiming(fileId,[0,i]) ax.plot(delta_t,dNdt,linestyle="steps-mid",label="gen=%.0f"%gen) i+=1 - else: - delta_t,dNdt = ReadTiming(fileId,[0,3]) - ax.plot(delta_t,dNdt,linestyle="steps-mid",label=fileId) ax.set_xscale('log') ax.set_yscale('log') @@ -81,12 +66,19 @@ def drawConvolveTiming(files): fig = figure() ax = fig.add_subplot(111) for fileId in files: - delta_t,dNdt1 = ReadTiming(fileId,[0,8]) # MeV band - delta_t,dNdt2 = ReadTiming(fileId,[0,9]) # GeV band - delta_t,dNdt3 = ReadTiming(fileId,[0,10]) # TeV band - ax.plot(delta_t,convolve(dNdt1,dNdt2,'same'),linestyle="steps-mid",label="MeV - GeV band") - ax.plot(delta_t,convolve(dNdt2,dNdt3,'same'),linestyle="steps-mid",label="GeV - TeV band") - ax.plot(delta_t,convolve(dNdt1,dNdt3,'same'),linestyle="steps-mid",label="MeV - TeV band") + delta_t,dNdt1 = ReadTiming(fileId,[0,2]) # MeV band + delta_t,dNdt2 = ReadTiming(fileId,[0,3]) # GeV band + delta_t,dNdt3 = ReadTiming(fileId,[0,4]) # TeV band + conv1 = convolve(dNdt1,dNdt2,'same') + conv2 = convolve(dNdt2,dNdt3,'same') + conv3 = convolve(dNdt1,dNdt3,'same') + ax.plot(delta_t,conv1/max(conv1),linestyle="steps-mid",label="MeV - GeV band") + ax.plot(delta_t,conv2/max(conv2),linestyle="steps-mid",label="GeV - TeV band") + ax.plot(delta_t,conv3/max(conv3),linestyle="steps-mid",label="MeV - TeV band") + print fileId, " =====================================" + print " MeV - GeV band: ",delta_t[r_[True, conv1[1:] > conv1[:-1]] & r_[conv1[:-1] > conv1[1:], True]] + print " GeV - TeV band: ",delta_t[r_[True, conv2[1:] > conv2[:-1]] & r_[conv2[:-1] > conv2[1:], True]] + print " MeV - TeV band: ",delta_t[r_[True, conv3[1:] > conv3[:-1]] & r_[conv3[:-1] > conv3[1:], True]] ax.set_xscale('log') ax.set_yscale('log') @@ -117,21 +109,26 @@ def drawDelay_vs_angle(fileId): ax = fig.add_subplot(111) nbBins = 100 - theta = ReadExtraFile(fileId,[4])*degre + generation, theta = ReadExtraFile(fileId,cols=[3,4]) + theta *= degre energy = ReadEnergy(fileId) delay = ReadTime(fileId) Emin = [1e-3,1e0,1e3] #GeV Emax = [1e0,1e3,1e5] #GeV - for n in arange(0,3,1): - cond= (energy>Emin[n]) & (energyEmin[n]) & (energy100) & (e<1e5) + e=e[cond] + g=g[cond] + f=f[cond] + p = ax11.plot(e,g,"-"+color[ind],label="z = "+lab) + ax11.plot(e,f,color=p[0].get_color(),linestyle='--') + + ind=ind+1 + +ax11.axvline(x=Ethreshold_gg(), ymin=0., ymax = 1., color='r', linewidth=2) + +ax11.legend(loc="best") +ax11.grid(b=True,which='major') +ax11.set_ylim([1e-4,1e4]) +ax11.set_xscale('log') +ax11.set_yscale('log') +ax11.set_xlabel("energy [GeV]") +ax11.set_ylabel("$\lambda_{\gamma\gamma}$ [Mpc]") + +plt.show() -- libgit2 0.21.2