diff --git a/Analysis.py b/Analysis.py index 2678453..b705273 100755 --- a/Analysis.py +++ b/Analysis.py @@ -140,7 +140,7 @@ for fileId in argv[1:]: 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)Elim) #& (abs(theta)Emin[n]) & (energy3: + fy=1 + i2-=1 + i1=i2 + else: + i1 = i2-1 + fy=(z-z_tab[i1])/(z_tab[i2]-z_tab[i1]) + j2 = searchsorted(E_tab,Egamma) + j1 = j2-1 + fx=(Egamma-E_tab[j1])/(E_tab[j2]-E_tab[j1]) + lambda_e = loadtxt("lambda_e.dat",unpack=True,usecols=[i1+1,i2+1,i1+6,i2+6]) + lambda11 = lambda_e[0,j1]*((1-fx)*(1-fy)) + lambda12 = lambda_e[1,j1]*((1-fx)*fy) + lambda21 = lambda_e[0,j2]*(fx*(1-fy)) + lambda22 = lambda_e[1,j2]*(fx*fy) + lambda_proper = lambda11 + lambda12 + lambda21 + lambda22 + lambda11 = lambda_e[2,j1]*((1-fx)*(1-fy)) + lambda12 = lambda_e[3,j1]*((1-fx)*fy) + lambda21 = lambda_e[2,j2]*(fx*(1-fy)) + lambda22 = lambda_e[3,j2]*(fx*fy) + lambda_comobile = lambda11 + lambda12 + lambda21 + lambda22 + return lambda_proper, lambda_comobile, 800.e3/Egamma #Mpc (from Durrer and Neronov 2013) diff --git a/Modules/Analytic.pyc b/Modules/Analytic.pyc index e9713a5..2bd40ab 100644 Binary files a/Modules/Analytic.pyc and b/Modules/Analytic.pyc differ diff --git a/Modules/Angle.py b/Modules/Angle.py index ea3c323..5d52a49 100644 --- a/Modules/Angle.py +++ b/Modules/Angle.py @@ -11,7 +11,6 @@ def angle_vs_energy(theta,energy,weight,nbBins=50): 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) @@ -23,9 +22,9 @@ def angle_vs_energy(theta,energy,weight,nbBins=50): def radial(theta,weight,nbBins=50): cond = theta!=0 - theta = 2*log10(theta[cond]) + theta = log10(theta[cond]) weight = weight[cond] - dn,angle2 = histogram(theta,nbBins,range=[2*log10(1e-3),2*log10(25)],weights=weight) + dn,angle2 = histogram(theta,nbBins,range=[log10(1e-3),log10(25)],weights=weight) angle2=10**angle2 #theta = theta**2 #dn,angle2 = histogram(theta,nbBins,range=[0,25],weights=weight) @@ -66,8 +65,8 @@ def drawRadial(files,all_source_spectrum=False,all_jets=False): ax.set_xscale('log') ax.set_yscale('log') ax.grid(b=True,which='major') - ax.set_xlabel("$\\theta^2$ [deg$^2$]") - ax.set_ylabel("$dN/d\\theta^2$ [arbitrary]") + ax.set_xlabel("$\\theta$ [deg]") + ax.set_ylabel("$dN/d\\theta$ [arbitrary]") show() @@ -97,7 +96,7 @@ def drawAngle_vs_energy(files,all_source_spectrum=False,all_jets=False,PlotAnaly if PlotAnalytic: yfit = Analytic_theta(ener,fileId) - ax1.plot(ener,yfit,'--',color=p[0].get_color(),linewidth=2,label="analytic") + ax1.plot(ener,yfit,'--',color=p[0].get_color(),linewidth=2) ax1.legend(loc="best") ax1.set_xscale('log') diff --git a/Modules/Angle.pyc b/Modules/Angle.pyc index 6ca61eb..9edcd34 100644 Binary files a/Modules/Angle.pyc and b/Modules/Angle.pyc differ diff --git a/Modules/Constants.pyc b/Modules/Constants.pyc index 7062e5a..3bd6752 100644 Binary files a/Modules/Constants.pyc and b/Modules/Constants.pyc differ diff --git a/Modules/Read.pyc b/Modules/Read.pyc index 5187303..e4d6b42 100644 Binary files a/Modules/Read.pyc and b/Modules/Read.pyc differ diff --git a/Modules/Spectrum.py b/Modules/Spectrum.py index 44ff8e1..ef42fd0 100644 --- a/Modules/Spectrum.py +++ b/Modules/Spectrum.py @@ -40,10 +40,11 @@ def drawSpectrum(files,all_source_spectrum=False,all_jets=False,PlotAnalytic=Fal Output: graph spectrum normalize to 1 source photon ''' fig = figure() - gs = gridspec.GridSpec(2, 1, height_ratios=[3,1]) - fig.subplots_adjust(hspace=0.001) - ax1 = fig.add_subplot(gs[0]) - ax2 = fig.add_subplot(gs[1],sharex=ax1) + #gs = gridspec.GridSpec(2, 1, height_ratios=[3,1]) + #fig.subplots_adjust(hspace=0.001) + #ax1 = fig.add_subplot(gs[0]) + #ax2 = fig.add_subplot(gs[1],sharex=ax1) + ax1 = fig.add_subplot(111) for fileId in files: if all_source_spectrum: @@ -57,9 +58,9 @@ def drawSpectrum(files,all_source_spectrum=False,all_jets=False,PlotAnalytic=Fal 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()) + #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: @@ -73,23 +74,24 @@ def drawSpectrum(files,all_source_spectrum=False,all_jets=False,PlotAnalytic=Fal 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()) + #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 else: # read files - Es,Fs = ReadSourceSpectrum(fileId,[0,3]) - p = ax1.plot(Es,Fs,linestyle=':') + #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,color=p[0].get_color(),drawstyle='steps-mid',label=fileId) + 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()) + #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 PlotAnalytic: minEner=min(ener) @@ -97,16 +99,18 @@ def drawSpectrum(files,all_source_spectrum=False,all_jets=False,PlotAnalytic=Fal 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) ax1.legend(loc="best") ax1.set_xscale('log') ax1.set_yscale('log') ax1.grid(b=True,which='major') ax1.set_ylabel("$E^2 dN/dE$ [GeV]") - ax2.set_xscale('log') - ax2.grid(b=True,which='major') - ax2.set_xlabel("energy [GeV]") - xticklabels = ax1.get_xticklabels() - setp(xticklabels, visible=False) + ax1.set_xlabel("energy [GeV]") + #ax2.set_xscale('log') + #ax2.grid(b=True,which='major') + #ax2.set_xlabel("energy [GeV]") + #xticklabels = ax1.get_xticklabels() + #setp(xticklabels, visible=False) show() diff --git a/Modules/Spectrum.pyc b/Modules/Spectrum.pyc index 3fbaaae..2593ae4 100644 Binary files a/Modules/Spectrum.pyc and b/Modules/Spectrum.pyc differ diff --git a/Modules/Timing.py b/Modules/Timing.py index eb84f1f..d34759b 100644 --- a/Modules/Timing.py +++ b/Modules/Timing.py @@ -1,4 +1,4 @@ -from numpy import histogram, log10, shape, arange, pi +from numpy import histogram, log10, shape, arange, pi, convolve from matplotlib.pyplot import figure, show from Read import ReadTime, ReadTiming, ReadExtraFile, ReadProfile, ReadEnergy from Constants import yr, degre @@ -11,39 +11,59 @@ def timing(time,weight,nbBins=100): Input (optional): generation (default = all) Output: time delay (sec), flux ''' - time=time[time>0] - weight=weight[time>0] + cond = time>0 + time = time[cond] + time=log10(time) + weight=weight[cond] - time = log10(time*yr) - - dN,dt=histogram(time,nbBins,weights=weight) + #dN,dt=histogram(time,nbBins,range=[-1,1e17],weights=weight) + dN,dt=histogram(time,nbBins,range=[-1,17],weights=weight) + dt = 10**dt timecenter=(dt[1:nbBins+1]+dt[0:nbBins])/2 binSize=dt[1:nbBins+1]-dt[0:nbBins] - dNdt=dN/(max(dN)*binSize) + dNdt=dN/binSize *timecenter - return 10**timecenter, dNdt + return timecenter, dNdt -def drawTiming(files,all_source_spectrum=False,all_jets=False): +def drawTiming(files,plot_case=0): ''' 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 all_source_spectrum: + 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 all_jets: + 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 + 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) + i+=1 + elif plot_case == 4: + i=11 + generation = ReadExtraFile(fileId,[3]) + for gen in list(set(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) @@ -53,10 +73,40 @@ def drawTiming(files,all_source_spectrum=False,all_jets=False): ax.grid(b=True,which='major') ax.legend(loc="best") ax.set_xlabel("Time delay [s]") - ax.set_ylabel("$t\ dN/dt$ [arbitrary units]") + ax.set_ylabel("$t.dN/dt$ [arbitrary units]") + + show() + +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") + + ax.set_xscale('log') + ax.set_yscale('log') + ax.grid(b=True,which='major') + ax.legend(loc="best") + ax.set_xlabel("Time delay [s]") + ax.set_ylabel("cross correlation [arbitrary units]") show() +def delay_vs_theta(theta,delay,nbBins=100): + cond= (theta!=0) + theta=log10(theta[cond]) + dt,dtheta=histogram(theta,nbBins,weights=delay[cond]) + dN,dtheta=histogram(theta,nbBins) + dtheta=10**dtheta + thetacenter=(dtheta[1:nbBins+1]+dtheta[0:nbBins])/2 + dt=dt/dN + return dt, thetacenter + def drawDelay_vs_angle(fileId): ''' Plot angle versus time delay, generation by generation @@ -67,38 +117,37 @@ def drawDelay_vs_angle(fileId): ax = fig.add_subplot(111) nbBins = 100 - theta = ReadExtraFile(fileId,[4]) + theta = ReadExtraFile(fileId,[4])*degre energy = ReadEnergy(fileId) delay = ReadTime(fileId) - for n in arange(0,8,2): - Emin = 10**(3-n) #GeV - Emax = 10**(5-n) #GeV - label = "%1.0e Gev"%Emin+"< $E_\\gamma$ < %1.0e Gev"%Emax - cond= (energy>Emin) & (energyEmin[n]) & (energy