From 484df0323a4ce10dffb6fa0eab80bfc54ed0dd08 Mon Sep 17 00:00:00 2001 From: Thomas Fitoussi Date: Fri, 4 Sep 2015 14:28:08 +0200 Subject: [PATCH] Correction on the theta^2 distribution and on the analytic expression of theta vs E --- Analysis.py | 28 +++++++++++++++------------- Modules/Analytic.py | 9 ++++----- Modules/Analytic.pyc | Bin 4508 -> 0 bytes Modules/Angle.py | 44 +++++++++++++++++++++++++++++--------------- Modules/Angle.pyc | Bin 3340 -> 0 bytes 5 files changed, 48 insertions(+), 33 deletions(-) diff --git a/Analysis.py b/Analysis.py index 99e27f8..2678453 100755 --- a/Analysis.py +++ b/Analysis.py @@ -6,7 +6,7 @@ from Modules.Read import ReadEnergy, ReadTime, ReadExtraFile, ReadProfile from Modules.Read import ReadMomentumAngle, ReadPositionAngle, resultDirectory from Modules.Spectrum import spectrum from Modules.Map import computeMap -from Modules.Angle import angle_vs_energy +from Modules.Angle import angle_vs_energy, radial from Modules.Timing import timing from Modules.Constants import degre @@ -29,8 +29,9 @@ for fileId in argv[1:]: energy = ReadEnergy(fileId) weightini, generation, theta_arrival, Esource, dir_source = ReadExtraFile(fileId,[2,3,4,5,6]) #weightini, generation, theta_arrival, Esource = ReadExtraFile(fileId,[2,3,4,5]) - nbPhotonsEmitted=ReadProfile(fileId,[3]) - weightini /= nbPhotonsEmitted + n_phot, n_lept = ReadProfile(fileId,[3,5]) + ratio = n_phot#/n_lept + weightini /= ratio theta_arrival*= degre thetaDir,phiDir = ReadMomentumAngle(fileId)*degre @@ -68,7 +69,7 @@ for fileId in argv[1:]: nbBins = 100 # draw source spectrum Es=array(list(set(Esource))) - Ws= (Es/min(Es))**(1-powerlaw_index) /nbPhotonsEmitted + Ws= (Es/min(Es))**(1-powerlaw_index) / ratio Es,Fs = spectrum(Es,Ws,nbBins=nbBins) Es=Es[:,newaxis] Fs=Fs[:,newaxis] @@ -92,8 +93,10 @@ for fileId in argv[1:]: Spectrum = append(Spectrum,flux_0,axis=1) # IMAGING, RADIAL DISTRIBUTION AND ANGLE VERSUS ENERGY ======================# - nbBins = 60 - theta2,dndtheta2,ener,angle = angle_vs_energy(theta_arrival,energy,weight,nbBins) + 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] @@ -128,7 +131,7 @@ for fileId in argv[1:]: Jet_Opening=[180,60,30] # degre (180 <=> isotrop) powerlaw_index = 2 Elim = 1e-1 # GeV - thetalim = 50 # degre + thetalim = 20 # degre # apply source spectrum weight_source = (Esource/min(Esource))**(1-powerlaw_index) @@ -160,7 +163,7 @@ for fileId in argv[1:]: nbBins = 100 # draw source spectrum Es=array(list(set(Esource))) - Ws= (Es/min(Es))**(1-powerlaw_index) /nbPhotonsEmitted + Ws= (Es/min(Es))**(1-powerlaw_index) / ratio Es,Fs = spectrum(Es,Ws,nbBins=nbBins) Es=Es[:,newaxis] Fs=Fs[:,newaxis] @@ -184,9 +187,10 @@ for fileId in argv[1:]: Spectrum = append(Spectrum,flux_0,axis=1) # IMAGING, RADIAL DISTRIBUTION AND ANGLE VERSUS ENERGY ======================# - nbBins = 60 + nbBins = 50 computeMap(theta,phi,weight,energy,fileId,jet_opening_angle,nbBins,borne=[thetalim,thetalim]) - theta2,dndtheta2,ener,angle = angle_vs_energy(theta_arrival,energy,weight,nbBins) + ener,angle = angle_vs_energy(theta_arrival,energy,weight,nbBins) + theta2,dndtheta2 = radial(theta_arrival,weight,nbBins) theta2=theta2[:,newaxis] dndtheta2=dndtheta2[:,newaxis] ener=ener[:,newaxis] @@ -201,9 +205,7 @@ for fileId in argv[1:]: Radial = append(Radial,dndtheta2,axis=1) Angle_Energy = append(Angle_Energy,angle,axis=1) - #=============================================================================# - # TIME DISTRIBUTION AND TIME DELAY VERSUS ANGLE - #=============================================================================# + # TIME DISTRIBUTION AND TIME DELAY VERSUS ANGLE =============================# nbBins = 100 delta_t,dNdt = timing(time,weight,nbBins) delta_t=delta_t[:,newaxis] diff --git a/Modules/Analytic.py b/Modules/Analytic.py index 4e73cf4..1bc6a73 100644 --- a/Modules/Analytic.py +++ b/Modules/Analytic.py @@ -16,9 +16,8 @@ def Analytic_theta(E_ic,fileId): delta_to_theta=lambda_gg(Esource)/Dsource RL0=RL(Esource/2,B) Dic0=Dic(Esource/2) - Ee2= E_ic /Eic(1) #GeV^2 - delta=Dic0/(2*RL0)*((Esource/2)**2 /Ee2 -1) - return abs(arcsin(delta_to_theta*sin(delta))) + delta=Dic0/(2*RL0)*((Esource/(2*Ee(E_ic)))**2 -1) + return abs(arcsin(delta_to_theta*sin(delta)))*degre def Analytic_delay_vs_theta(theta,fileId): Esource,distSource=ReadProfile(fileId,[0,2]) #GeV #Mpc @@ -60,7 +59,7 @@ def Eic(Ee=Ee_default): return 4*Ecmb*Ee**2/(3*me**2)*1e3 #GeV def Ee(Egamma=Egamma_default): - return sqrt((3*me**2)/(4*Ecmb)*1e-3 *Egamma) #GeV + return me*sqrt((3*Egamma*1e-3 )/(4*Ecmb)) #GeV def tIC(): return lambdaIC()/(c*yr/Mpc) #yr @@ -87,4 +86,4 @@ def Ethreshold_gg(): return (me)**2/Eebl *1e-3 #GeV def lambda_gg(Egamma=1): # Egamma (GeV) - return 800. /(Egamma)#Mpc (from Durrer and Neronov 2013) + return 800. /(Egamma) *1e3#Mpc (from Durrer and Neronov 2013) diff --git a/Modules/Analytic.pyc b/Modules/Analytic.pyc index bcbb67b..e9713a5 100644 Binary files a/Modules/Analytic.pyc and b/Modules/Analytic.pyc differ diff --git a/Modules/Angle.py b/Modules/Angle.py index 6ab6cfe..ea3c323 100644 --- a/Modules/Angle.py +++ b/Modules/Angle.py @@ -19,20 +19,22 @@ def angle_vs_energy(theta,energy,weight,nbBins=50): enercenter = (ener[1:nbBins+1]+ener[0:nbBins])/2 angle /= nb - # RADIAL DISTRIBUTION ========================================================# + return enercenter, angle + +def radial(theta,weight,nbBins=50): cond = theta!=0 - #theta = 2*log10(theta[cond]) - theta = theta[cond]**2 + theta = 2*log10(theta[cond]) weight = weight[cond] - #dn,angle2 = histogram(theta,nbBins,range=[2*log10(1e-6),2*log10(90)],weights=weight) - dn,angle2 = histogram(theta,nbBins,range=[0,20],weights=weight) - #angle2=10**angle2 + dn,angle2 = histogram(theta,nbBins,range=[2*log10(1e-3),2*log10(25)],weights=weight) + angle2=10**angle2 + #theta = theta**2 + #dn,angle2 = histogram(theta,nbBins,range=[0,25],weights=weight) thetacenter = (angle2[1:nbBins+1]+angle2[0:nbBins])/2 dtheta = angle2[1:nbBins+1]-angle2[0:nbBins] dndtheta = dn/dtheta dndtheta /= max(dndtheta) - return thetacenter, dndtheta, enercenter, angle + return thetacenter, dndtheta def drawRadial(files,all_source_spectrum=False,all_jets=False): ''' @@ -61,15 +63,15 @@ def drawRadial(files,all_source_spectrum=False,all_jets=False): ax.plot(theta,dndtheta2,drawstyle='steps-mid',label=fileId) ax.legend(loc="best") - #ax.set_xscale('log') + 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$ [deg$^{-1}$]") + ax.set_ylabel("$dN/d\\theta^2$ [arbitrary]") show() -def drawAngle_vs_energy(files): +def drawAngle_vs_energy(files,all_source_spectrum=False,all_jets=False,PlotAnalytic=False): ''' Plot arrival angle versus energy + analytic expression Input: list of directories @@ -79,18 +81,30 @@ def drawAngle_vs_energy(files): ax1 = fig.add_subplot(111) for fileId in files: - ener,angle = ReadAngleVsEnergy(fileId,[0,1]) - p=ax1.plot(ener,angle,drawstyle='steps-mid',label=fileId) + 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,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) - yfit = Analytic_theta(ener,fileId) - ax1.plot(ener,yfit,'--',color="m",linewidth=2,label="analytic") + if PlotAnalytic: + yfit = Analytic_theta(ener,fileId) + ax1.plot(ener,yfit,'--',color=p[0].get_color(),linewidth=2,label="analytic") ax1.legend(loc="best") ax1.set_xscale('log') ax1.set_yscale('log') ax1.set_ylim([0,180]) ax1.grid(b=True,which='major') - ax1.set_ylabel("$\\theta$ [rad]") + ax1.set_ylabel("$\\theta$ [deg]") ax1.set_xlabel("energy [GeV]") show() diff --git a/Modules/Angle.pyc b/Modules/Angle.pyc index 125c20d..6ca61eb 100644 Binary files a/Modules/Angle.pyc and b/Modules/Angle.pyc differ -- libgit2 0.21.2