diff --git a/Analysis.py b/Analysis.py new file mode 100755 index 0000000..bfd9a53 --- /dev/null +++ b/Analysis.py @@ -0,0 +1,37 @@ +#!/bin/python + +from sys import argv +from numpy import shape +from Modules.Spectrum import drawSpectrum +from Modules.Angle import drawAngle, drawRadial +from Modules.Timing import drawTiming +from Modules.Generation import drawGeneration +from Modules.Map import drawMap, drawMapHisto + +if shape(argv)[0] < 3: + print "Give at least 2 arguments" + exit() + +if argv[1] == "spectrum": + drawSpectrum(argv[2:]) + +elif argv[1] == "angle": + drawAngle(argv[2:]) + +elif argv[1] == "radial": + drawRadial(argv[2:]) + +elif argv[1] == "timing": + drawTiming(argv[2:]) + +elif argv[1] == "generation": + drawGeneration(argv[2:]) + +elif argv[1] == "maphisto": + drawMapHisto(argv[2:]) + +elif argv[1] == "map": + drawMap(argv[2:]) + +else: + print "bad 1st argument" diff --git a/Angle.py b/Angle.py deleted file mode 100755 index 689e8d7..0000000 --- a/Angle.py +++ /dev/null @@ -1,135 +0,0 @@ -#!/usr/bin/python - -from sys import argv -import matplotlib.pyplot as plt -from matplotlib import gridspec -from Constantes import * - -if argv[1] == "EGMF": - fileName="Results_EGMF" -elif argv[1] == "EBL": - fileName="Results_EBL" - B=10**(-15) - Model=["Kneiske and Doll - 'best fit'","Kneiske and Doll - 'lower limit'", - "Fraceschini","Finke and Al","Gilmore and Al","Dominguez and Al"] -else: - print "Used: ./Angle_distribution.py arg1 ind1 ind2 ..." - print " arg1 = EGMF or EBL " - print " ind1 = file number \n" - print "Examples: " - print " to study EGMF14: ./Angle_distribution.py EGMF 14 " - print " to study EBL1 and EBL2: ./Angle_distribution.py EBL 1 2 \n" - exit() - -nbBins = 100 -Eliminf = 1e0 #GeV -Elimsup = 1e5 #GeV - -Esource=100 #TeV -Dsource=135 #Mpc -delta_to_theta=(lambda_gg/Esource*1e3)/Dsource - -def fit_theta(E,B): - RL0=RL*(Esource/2)*(1e-17/B) - Dic0=Dic/(Esource/2) - delta=Dic0/(2*RL0)*((Esource/2)**2 *Eic/E-1) - return abs(arcsin(delta_to_theta*sin(delta))*degre) - -fig2 = plt.figure() -ax2 = fig2.add_subplot(111) - -fig1 = plt.figure() -gs = gridspec.GridSpec(2, 1, height_ratios=[4,1]) -fig1.subplots_adjust(hspace=0.001) -ax11 = fig1.add_subplot(gs[0]) -ax12 = fig1.add_subplot(gs[1],sharex=ax11) - -color=['b','r','g','c','m','y'] - -ind = 0 -for fileId in argv[2:]: - energy,dirx,diry,dirz = loadtxt(fileName+fileId+"/Results_momentum",unpack=True,usecols=[0,1,2,3]) - posx,posy,posz = loadtxt(fileName+fileId+"/Results_position",unpack=True,usecols=[1,2,3]) - charge,weight,gen=loadtxt(fileName+fileId+"/Results_extra",unpack=True,usecols=[0,2,3]) - - hyp=sqrt(posx**2+posy**2+posz**2) - posx=posx/hyp - posy=posy/hyp - posz=posz/hyp - hyp=sqrt(dirx**2+diry**2+dirz**2) - dirx=dirx/hyp - diry=diry/hyp - dirz=dirz/hyp - costheta=dirx*posx+diry*posy+dirz*posz - cond=(costheta<=1)&(costheta>=-1)&(energy>Eliminf)&(energy> lambda_B -Delta2=sqrt(Dic*lambda_B)/RL*degre #degre if Dic << lambda_B - -# Threshold energy when Dic = RL -Ethreshold_ic = Eic*Dic/RL # GeV -# Pair production -Ethreshold_gg=(me)**2/Eebl *1e-3 #GeV -lambda_gg=0.8 # (E_gamma/1TeV)^-1 Gpc (from Durrer and Neronov 2013) - -# usefull integrand -def comobileTime(z): - return -1/(H0*(1+z)*sqrt(omegaM*(1+z)**3+omegaK*(1+z)**2+omegaL)) - -def distPhoton(z): - return -c/(H0*a0*sqrt(omegaM*(1+z)**3+omegaK*(1+z)**2+omegaL)) - -def distLepton(z,E): - beta = sqrt(1-m**2*c**4/E**2) - return -beta*c/(H0*a0*sqrt(omegaM*(1+z)**3+omegaK*(1+z)**2+omegaL)) diff --git a/Constantes.pyc b/Constantes.pyc deleted file mode 100644 index 82238ce..0000000 Binary files a/Constantes.pyc and /dev/null differ diff --git a/Generation.py b/Generation.py deleted file mode 100755 index 271243e..0000000 --- a/Generation.py +++ /dev/null @@ -1,34 +0,0 @@ -#!/usr/bin/python - -from sys import argv -import numpy as np -import matplotlib.pyplot as plt - -if len(argv)<2: - print "Please give at least 1 file id \n For example: 14 for Results.EGMF14 \n or 14 15 for EGMF14 and EGMF15" - exit() - -# figure: energy distribution -#============================= -nbBins = 7 -color=['b','g','r','c','m','y'] - -fig = plt.figure() -ax = fig.add_subplot(111) - -ind = 0 -for fileId in argv[1:]: - weight,nbGen=np.loadtxt("Results_EGMF"+fileId+"/Results_extra",unpack=True,usecols=[2,3]) - - plt.hist(nbGen,nbBins,range=[2,9],log=1,facecolor=color[ind],alpha=.5, weights=weight, - label="$10^{-%.0f"%float(fileId)+"}$Gauss") - ind=ind+1 - -xtick=np.arange(9) -ax.legend(loc='best') -ax.set_xticks(xtick+0.5) -ax.set_xticklabels(xtick) -ax.set_xlabel("Generation") -ax.set_ylabel("Number of events") - -plt.show() diff --git a/Map_arrival.py b/Map_arrival.py deleted file mode 100755 index e489e9f..0000000 --- a/Map_arrival.py +++ /dev/null @@ -1,83 +0,0 @@ -#!/usr/bin/python - -from sys import argv -from numpy import * -import matplotlib.pyplot as plt -from matplotlib.colors import LogNorm -from mpl_toolkits.mplot3d.axes3d import Axes3D - -if len(argv)<2: - print "Please give at least 1 file id \n For example: 14 for Results.EGMF14 \n or 14 15 for EGMF14 and EGMF15" - exit() - -# figure: position of the arriving particles in space -> size circle = energy -#============================================================================= -nbBins = 100 -degre= 180/pi -borne=5 # degree -Elim=1e-1 # GeV - -fig = plt.figure() - -nbplots = len(argv)-1 -ind = 1 -for fileId in argv[1:]: - energy,thetaDir,phiDir = loadtxt("Results_EGMF"+fileId+"/Results_momentum",unpack=True,usecols=[0,4,5]) - thetaPos,phiPos = loadtxt("Results_EGMF"+fileId+"/Results_position",unpack=True,usecols=[4,5]) - charge,weight= loadtxt("Results_EGMF"+fileId+"/Results_extra",unpack=True,usecols=[0,2]) - - thetaPos = thetaPos*degre - phiPos = phiPos*degre - thetaDir = thetaDir*degre - phiDir = phiDir*degre - theta = thetaDir - thetaPos - phi = phiDir -phiPos - - cond=(energy>Elim) & (abs(theta)", shape(Energy)[0], "particles" - plotid = 100+nbplots*10+ind - - #Emax = max(Energy) - #Emin = min(Energy) - #area = pi * (15 * Energy/Emax)**2 +15 - #colors = Energy/Emax - - #ax = fig.add_subplot(plotid) - #sax = ax.scatter(theta, phi, s=area, c=colors, alpha=0.3, marker='o') - #ax.set_xlim((-borne,borne)) - #ax.set_ylim((-borne,borne)) - #ax.set_xlabel("$\\theta$ in degree") - #ax.set_ylabel("$\\phi$ in degree") - #ax.set_title("EGMF: $10^{-%.0f"%float(fileId)+"}$Gauss") - - #if Emin < 1e-3: - # liminf = "%.0f"%(Emin*1e6)+"keV" - #elif Emin < 1: - # liminf = "%.0f"%(Emin*1e3)+"MeV" - #else: - # liminf = "%.0f"%(Emin)+"GeV" - #limsup = "%.0f"%(Emax*1e-3)+"TeV" - #cbar = fig.colorbar(sax, ticks=[0.01,1]) - #cbar.ax.set_yticklabels([liminf,limsup]) - - ########################################### - - ax = fig.add_subplot(plotid) - H,xedges,yedges,im = ax.hist2d(theta,phi,nbBins,weights=weight/shape(weight),norm=LogNorm()) - extent = [yedges[0], yedges[-1], xedges[0], xedges[-1]] - cbar=fig.colorbar(im, ax=ax) - cbar.ax.set_ylabel("counts") - ax.set_xlim((-borne,borne)) - ax.set_ylim((-borne,borne)) - ax.set_xlabel("$\\theta$ in degree") - ax.set_ylabel("$\\phi$ in degree") - ax.set_title("EGMF: $10^{-%.0f"%float(fileId)+"}$Gauss") - - ind=ind+1 - -plt.show() diff --git a/Modules/Analytic.py b/Modules/Analytic.py new file mode 100644 index 0000000..4f88680 --- /dev/null +++ b/Modules/Analytic.py @@ -0,0 +1,15 @@ +from Constants import me, Ecmb, Eic, lambda_gg, degre, RL, Dic +from numpy import sqrt, abs, sin, arcsin + +def Analytic_flux(E_ic): + return me*1e3/4*sqrt(3e9/Ecmb)*1e-9*sqrt(E_ic) #GeV + +def Analytic_theta(E,fileId): + Esource=100 #TeV + Dsource=135 #Mpc + B=1e-15 #Gauss + delta_to_theta=(lambda_gg/Esource*1e3)/Dsource + RL0=RL*(Esource/2)*(1e-17/B) + Dic0=Dic/(Esource/2) + delta=Dic0/(2*RL0)*((Esource/2)**2 *Eic/E-1) + return abs(arcsin(delta_to_theta*sin(delta))*degre) diff --git a/Modules/Analytic.pyc b/Modules/Analytic.pyc new file mode 100644 index 0000000..ddcedc9 Binary files /dev/null and b/Modules/Analytic.pyc differ diff --git a/Modules/Angle.py b/Modules/Angle.py new file mode 100644 index 0000000..f4799bd --- /dev/null +++ b/Modules/Angle.py @@ -0,0 +1,103 @@ +from numpy import sqrt, arccos, shape, log10, histogram +from matplotlib.pyplot import figure, show, ylim, setp +from matplotlib import gridspec +from Read import ReadMomentum, ReadPosition, ReadEnergy, ReadWeight +from Analytic import Analytic_theta +from Constants import degre + +Eliminf = 1e0 #GeV +Elimsup = 1e5 #GeV + +def compute_theta(fileName): + position = ReadPosition(fileName) + momentum = ReadMomentum(fileName) + energy = ReadEnergy(fileName) + weight = ReadWeight(fileName) + + hyp=sqrt(position[0]**2+position[1]**2+position[2]**2) + position /= hyp + hyp=sqrt(momentum[0]**2+momentum[1]**2+momentum[2]**2) + momentum /= hyp + + costheta=position[0]*momentum[0]+position[1]*momentum[1]+position[2]*momentum[2] + cond=(costheta<=1)&(costheta>=-1)&(energy>Eliminf)&(energy", shape(energy)[0], "events" + + return energy, weight, theta + +def angle_vs_energy(fileName): + energy,weight,theta = compute_theta(fileName) + nbBins = 40 + + energy =log10(energy) + angle,ener =histogram(energy,nbBins,weights=theta) + nb,ener =histogram(energy,nbBins) + ener = 10**ener + enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2 + angle = angle/nb + + return enercenter, angle + +def drawAngle(files): + fig = figure() + gs = gridspec.GridSpec(2, 1, height_ratios=[4,1]) + fig.subplots_adjust(hspace=0.001) + ax1 = fig.add_subplot(gs[0]) + ax2 = fig.add_subplot(gs[1],sharex=ax1) + + for fileId in files: + energy,angle = angle_vs_energy(fileId) + p=ax1.plot(energy,angle,'.',label=fileId) + + yfit = Analytic_theta(energy,fileId) + ax1.plot(energy,yfit,'--',color=p[0].get_color(),linewidth=2,label="analytic") + + error=angle/yfit-1 + ax2.plot(energy,error,'+',color=p[0].get_color()) + + 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$ [degre]") + ax2.set_xscale('log') + ax2.grid(b=True,which='major') + ax2.set_xlabel("energy [GeV]") + ax2.set_ylabel("relative error") + step=0.2*(max(error)-min(error)) + ylim(min(error)-step, max(error)+step) + xticklabels = ax1.get_xticklabels() + setp(xticklabels, visible=False) + + show() + +def radial(fileName): + nbBins = 20 + energy,weight,theta = compute_theta(fileName) + + dn,angle=histogram(theta**2,nbBins,range=[0,0.25],weights=weight) + thetacenter=(angle[1:nbBins+1]+angle[0:nbBins])/2 + dtheta=angle[1:nbBins+1]-angle[0:nbBins] + dndtheta2=dn/dtheta + + return thetacenter, dndtheta2 + +def drawRadial(files): + fig = figure() + ax = fig.add_subplot(111) + + for fileId in files: + theta,dndtheta2=radial(fileId) + ax.plot(theta,dndtheta2,drawstyle='steps-mid',label=fileId) + + ax.legend(loc="best",title="%3d - %3d GeV"%(Eliminf,Elimsup)) + ax.set_yscale('log') + ax.grid(b=True,which='major') + ax.set_xlabel("$\\theta^2$ [deg$^2$]") + ax.set_ylabel("$dN/d\\theta^2$ [deg$^{-2}$]") + + show() diff --git a/Modules/Angle.pyc b/Modules/Angle.pyc new file mode 100644 index 0000000..873caf4 Binary files /dev/null and b/Modules/Angle.pyc differ diff --git a/Modules/Constants.py b/Modules/Constants.py new file mode 100644 index 0000000..61741d9 --- /dev/null +++ b/Modules/Constants.py @@ -0,0 +1,69 @@ +#!/bin/python + +from numpy import pi, sqrt +from scipy.special import zeta + +e=4.8032068e-10 # esu (cgs units) +qe=1.602176565e-12 # erg +m=9.1093897e-28 # g +me=510.998928 # kev +c=2.99792458e10 # cm.s-1 +k=1.380658E-16 # erg.K-1 +h=6.62606957e-27 # erg.s +hb=1.05457266E-27 # erg.s +alpha=1./137.035999679 +lambdaC=hb/(m*c) +Mpc=(3.0856776e+16)*1e8 # Mpc to cm +erg_to_GeV=1/(1.602176565e-3) +degre=180/pi +yr=3600*24*365.25 # s + +# Cosmology +a0=1. +H0=67.8*1e5/(Mpc) # s-1 +omegaM = 0.3 +omegaK = 0 +omegaL = 0.7 + +# CMB +Tcmb=2.7255 # K +Tcmbp=(k*Tcmb)/(m*c*c) # Adimentionnal thermal energy +nCMB=(16*pi*zeta(3,1))*(k*Tcmb/(h*c))**3 # cm^-3 +Ecmb=(pi**4/(30*zeta(3,1)))*Tcmbp*me*1e3 #eV +rhoCMB=Ecmb*nCMB #eV.cm^-3 + +# EBL +Eebl=1 #eV + +# Particles physics +r0=e*e/(m*c*c) +sigmaT=8.*pi*(r0**2)/3. # Thomson cross section + +# Compton accumulation +Compton_threshold = 0.05 +ECompton_threshold = Compton_threshold/(4/3*2.7*Tcmbp) *me*1e-6 #GeV + +# Value for analytic expression +Ee=1e3 #GeV +B=1e-17 #Gauss +lambda_B=0.3 #Mpc + +# Compton scattering +Dic= 3*(me*1e-6)**2/(4*sigmaT*rhoCMB*1e-9*Ee) /Mpc #Mpc +lambdaIC=1/(nCMB*sigmaT*Mpc) #Mpc +Eic= 4*Ecmb*Ee**2/(3*me**2)*1e3 #GeV +tIC=lambdaIC/(c*yr/Mpc) #yr + +# Larmor radius +RL=(Ee/erg_to_GeV)/(e*B) /Mpc #Mpc +# Magnetic deflection +delta=lambdaIC/RL*degre +Delta1=Dic/RL*degre #degre if Dic >> lambda_B +Delta2=sqrt(Dic*lambda_B)/RL*degre #degre if Dic << lambda_B + +# Threshold energy when Dic = RL +Ethreshold_ic = Eic*Dic/RL # GeV +# Pair production +Ethreshold_gg=(me)**2/Eebl *1e-3 #GeV +lambda_gg=0.8 # (E_gamma/1TeV)^-1 Gpc (from Durrer and Neronov 2013) + diff --git a/Modules/Constants.pyc b/Modules/Constants.pyc new file mode 100644 index 0000000..4e30c47 Binary files /dev/null and b/Modules/Constants.pyc differ diff --git a/Modules/Generation.py b/Modules/Generation.py new file mode 100644 index 0000000..b038938 --- /dev/null +++ b/Modules/Generation.py @@ -0,0 +1,23 @@ +from numpy import shape, arange +from matplotlib.pyplot import figure, show, hist +from Read import ReadGeneration, ReadWeight + + +def drawGeneration(files): + fig = figure() + ax = fig.add_subplot(111) + nbBins = 7 + for fileId in files: + weight = ReadWeight(fileId) + generation = ReadGeneration(fileId) + print "file", fileId, "->", shape(generation)[0], "events" + hist(generation,nbBins,range=[2,9],log=1,alpha=.5, weights=weight,label=fileId) + + xtick=arange(9) + ax.legend(loc='best') + ax.set_xticks(xtick+0.5) + ax.set_xticklabels(xtick) + ax.set_xlabel("Generation") + ax.set_ylabel("Number of events") + + show() diff --git a/Modules/Generation.pyc b/Modules/Generation.pyc new file mode 100644 index 0000000..79de4ef Binary files /dev/null and b/Modules/Generation.pyc differ diff --git a/Modules/Integrand.py b/Modules/Integrand.py new file mode 100644 index 0000000..0b68ce2 --- /dev/null +++ b/Modules/Integrand.py @@ -0,0 +1,12 @@ +from Constants import H0, omegaK, omegaM, omegaL, c, m , a0 +from numpy import sqrt + +def comobileTime(z): + return -1/(H0*(1+z)*sqrt(omegaM*(1+z)**3+omegaK*(1+z)**2+omegaL)) + +def distPhoton(z): + return -c/(H0*a0*sqrt(omegaM*(1+z)**3+omegaK*(1+z)**2+omegaL)) + +def distLepton(z,E): + beta = sqrt(1-m**2*c**4/E**2) + return -beta*c/(H0*a0*sqrt(omegaM*(1+z)**3+omegaK*(1+z)**2+omegaL)) diff --git a/Modules/Integrand.pyc b/Modules/Integrand.pyc new file mode 100644 index 0000000..b31c093 Binary files /dev/null and b/Modules/Integrand.pyc differ diff --git a/Modules/Map.py b/Modules/Map.py new file mode 100644 index 0000000..7ae50d7 --- /dev/null +++ b/Modules/Map.py @@ -0,0 +1,88 @@ +from numpy import shape, pi +from matplotlib.pyplot import figure, show +from matplotlib.colors import LogNorm +from Read import ReadEnergy, ReadWeight, ReadMomentumAngle, ReadPositionAngle +from Constants import degre + +Elim=1e-1 # GeV +borne=5 # degree +nbBins = 100 + +def thetaphi(fileId): + energy = ReadEnergy(fileId) + weight = ReadWeight(fileId) + thetaDir,phiDir = ReadMomentumAngle(fileId)*degre + thetaPos,phiPos = ReadPositionAngle(fileId)*degre + + theta = thetaDir - thetaPos + phi = phiDir -phiPos + + cond=(energy>Elim) & (abs(theta)", shape(energy)[0], "events" + return energy, theta, phi, weight + +def drawMap(files): + print "Cut at", Elim, "GeV" + fig = figure() + nbplots = len(files) + + ind=1 + for fileId in files: + energy, theta, phi, weight = thetaphi(fileId) + + ax = fig.add_subplot(100+nbplots*10+ind) + + Emax = max(energy) + Emin = min(energy) + colors = energy/Emax + area = pi * (15 * colors)**2 +15 + sax = ax.scatter(theta, phi, c=colors, s=area, alpha=0.3, marker='o') + + if Emin < 1e-3: + liminf = "%.0f"%(Emin*1e6)+"keV" + elif Emin < 1: + liminf = "%.0f"%(Emin*1e3)+"MeV" + else: + liminf = "%.0f"%(Emin)+"GeV" + limsup = "%.0f"%(Emax*1e-3)+"TeV" + cbar = fig.colorbar(sax, ticks=[0.01,1]) + cbar.ax.set_yticklabels([liminf,limsup]) + + ax.set_xlim((-borne,borne)) + ax.set_ylim((-borne,borne)) + ax.set_xlabel("$\\theta$ [deg]") + ax.set_ylabel("$\\phi$ [deg]") + ax.set_title(fileId) + + ind=ind+1 + + show() + +def drawMapHisto(files): + print "Cut at", Elim, "GeV" + fig = figure() + nbplots = len(files) + + ind=1 + for fileId in files: + energy, theta, phi, weight = thetaphi(fileId) + + ax = fig.add_subplot(100+nbplots*10+ind) + H,xedges,yedges,im = ax.hist2d(theta,phi,nbBins,weights=weight/shape(weight),norm=LogNorm()) + extent = [yedges[0], yedges[-1], xedges[0], xedges[-1]] + cbar=fig.colorbar(im, ax=ax) + cbar.ax.set_ylabel("counts") + ax.set_xlim((-borne,borne)) + ax.set_ylim((-borne,borne)) + ax.set_xlabel("$\\theta$ [deg]") + ax.set_ylabel("$\\phi$ [deg]") + ax.set_title(fileId) + + ind=ind+1 + + show() diff --git a/Modules/Map.pyc b/Modules/Map.pyc new file mode 100644 index 0000000..97d5f08 Binary files /dev/null and b/Modules/Map.pyc differ diff --git a/Modules/Read.py b/Modules/Read.py new file mode 100644 index 0000000..b0e3976 --- /dev/null +++ b/Modules/Read.py @@ -0,0 +1,44 @@ +from numpy import loadtxt + +resultDirectory="Results/" +positionFile="/position" +momentumFile="/momentum" +extraFile="/extra" +ElmagFile="/spec_diff_test" + +def ReadTime(fileName): + return loadtxt(resultDirectory+fileName+positionFile,unpack=True,usecols=[0]) + +def ReadPosition(fileName): + return loadtxt(resultDirectory+fileName+positionFile,unpack=True,usecols=[1,2,3]) + +def ReadPositionAngle(fileName): + return loadtxt(resultDirectory+fileName+positionFile,unpack=True,usecols=[4,5]) + +def ReadEnergy(fileName): + return loadtxt(resultDirectory+fileName+momentumFile,unpack=True,usecols=[0]) + +def ReadMomentum(fileName): + return loadtxt(resultDirectory+fileName+momentumFile,unpack=True,usecols=[1,2,3]) + +def ReadMomentumAngle(fileName): + return loadtxt(resultDirectory+fileName+momentumFile,unpack=True,usecols=[4,5]) + +def ReadCharge(fileName): + return loadtxt(resultDirectory+fileName+extraFile,unpack=True,usecols=[0]) + +def ReadRedshift(fileName): + return loadtxt(resultDirectory+fileName+extraFile,unpack=True,usecols=[1]) + +def ReadWeight(fileName): + return loadtxt(resultDirectory+fileName+extraFile,unpack=True,usecols=[2]) + +def ReadGeneration(fileName): + return loadtxt(resultDirectory+fileName+extraFile,unpack=True,usecols=[3]) + +def ReadElmag(fileName): + Elmag=resultDirectory+fileName+ElmagFile + energy,fluxGamma,fluxLepton = loadtxt(Elmag,unpack=True,usecols=[0,1,2]) + energy *= 10**(-9) + flux=fluxGamma+fluxLepton + return energy,flux diff --git a/Modules/Read.pyc b/Modules/Read.pyc new file mode 100644 index 0000000..df7ef64 Binary files /dev/null and b/Modules/Read.pyc differ diff --git a/Modules/Spectrum.py b/Modules/Spectrum.py new file mode 100644 index 0000000..db2e3e0 --- /dev/null +++ b/Modules/Spectrum.py @@ -0,0 +1,45 @@ +from numpy import histogram, log10, shape, logspace +from matplotlib.pyplot import figure, show +from os.path import isfile +from Read import ReadEnergy, ReadWeight, ReadElmag, resultDirectory, ElmagFile +from Analytic import Analytic_flux + +def spectrum(fileName): + energy = ReadEnergy(fileName) + weight = ReadWeight(fileName) + print "file", fileName, "->", shape(energy)[0], "events" + nbBins=100 + energy =log10(energy) + flux,ener=histogram(energy,nbBins,weights=weight) + ener = 10**ener + 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 + +def drawSpectrum(files): + fig = figure() + ax = fig.add_subplot(111) + for fileId in files: + energy,flux = spectrum(fileId) + maxflux = 10 #shape(energy)[0] + flux /= maxflux + p=ax.plot(energy,flux,drawstyle='steps-mid',label=fileId) + + Elmag=resultDirectory+fileId+ElmagFile + if isfile(Elmag): + energy, flux = ReadElmag(fileId) + ax.plot(energy,flux,"-",color=p[0].get_color(),label="Elmag - "+fileId) + + xfit = logspace(log10(10**(-3)),log10(2*10**4),200) + yfit = Analytic_flux(xfit) + ax.plot(xfit,yfit,'--m',linewidth=2,label="Analytic ($\\sqrt{E}$)") + + ax.set_xscale('log') + ax.set_yscale('log') + ax.grid(b=True,which='major') + ax.legend(loc="best") + ax.set_xlabel("energy [GeV]") + ax.set_ylabel("$E^2 \\frac{dN}{dE}$ [GeV]") + + show() diff --git a/Modules/Spectrum.pyc b/Modules/Spectrum.pyc new file mode 100644 index 0000000..6ee2338 Binary files /dev/null and b/Modules/Spectrum.pyc differ diff --git a/Modules/Timing.py b/Modules/Timing.py new file mode 100644 index 0000000..e8be8fb --- /dev/null +++ b/Modules/Timing.py @@ -0,0 +1,43 @@ +from numpy import histogram, log10, shape +from matplotlib.pyplot import figure, show +from scipy.integrate import quad +from Read import ReadTime, ReadWeight +from Constants import yr +from Integrand import comobileTime + +def timing(fileName): + time=ReadTime(fileName) + weight=ReadWeight(fileName) + prop = float(shape(time[time<0])[0])/shape(time)[0] + print "file", fileName, "->", shape(time)[0], "events", "negative time:", shape(time[time<0]), "~", prop + + time=time[time>0] + weight=weight[time>0] + + nbBins = 100 + zSource = 0.0308 + tlim = quad(comobileTime,zSource,0)[0]/yr + time=log10(time/tlim) + + dN,dt=histogram(time,nbBins,weights=weight) + timecenter=(dt[1:nbBins+1]+dt[0:nbBins])/2 + binSize=dt[1:nbBins+1]-dt[0:nbBins] + Nmax=shape(time)[0] + dNdt=dN/(Nmax*binSize) + + return timecenter, dNdt + +def drawTiming(files): + fig = figure() + ax = fig.add_subplot(111) + for fileId in files: + time,dNdt = timing(fileId) + ax.plot(time[dNdt!=0],dNdt[dNdt!=0],linestyle="steps-mid",label=fileId) + + ax.set_yscale('log') + ax.grid(b=True,which='major') + ax.legend(loc="best") + ax.set_xlabel("Time [$\Delta t / t$ log scale]") + ax.set_ylabel("$t\ dN/dt$") + + show() diff --git a/Modules/Timing.pyc b/Modules/Timing.pyc new file mode 100644 index 0000000..2c0cdb3 Binary files /dev/null and b/Modules/Timing.pyc differ diff --git a/Modules/__init__.py b/Modules/__init__.py new file mode 100644 index 0000000..e69de29 --- /dev/null +++ b/Modules/__init__.py diff --git a/Modules/__init__.pyc b/Modules/__init__.pyc new file mode 100644 index 0000000..d592137 Binary files /dev/null and b/Modules/__init__.pyc differ diff --git a/Spectrum.py b/Spectrum.py deleted file mode 100755 index 9e086a6..0000000 --- a/Spectrum.py +++ /dev/null @@ -1,87 +0,0 @@ -#!/usr/bin/python - -from sys import argv -import matplotlib.pyplot as plt -from Constantes import * - -if argv[1] == "EGMF": - fileName="Results_EGMF" -elif argv[1] == "EBL": - fileName="Results_EBL" - B=10**(-15) - Model=["Kneiske and Doll - 'best fit'","Kneiske and Doll - 'lower limit'", - "Franceschini","Finke and Al","Gilmore and Al","Dominguez and Al"] -else: - print "Used: ./Energy_distribution.py arg1 ind1 ind2 ..." - print " arg1 = EGMF or EBL " - print " ind1 = file number \n" - print "Examples: " - print " to study EGMF14: ./Energy_distribution.py EGMF 14 " - print " to study EBL1 and EBL2: ./Energy_distribution.py EBL 1 2 \n" - exit() - -# figure: energy distribution -#============================= -nbBins = 100 -convert= 180/pi -color=['b','r','g','c','m','y'] - -fig = plt.figure() -ax = fig.add_subplot(111) - -energy= loadtxt(fileName+argv[2]+"/Results_momentum",unpack=True,usecols=[0]) -Edata =log10(energy) -Emax = max(Edata) -Emin = min(Edata) - -ind = 0 -for fileId in argv[2:]: - energy= loadtxt(fileName+fileId+"/Results_momentum",unpack=True,usecols=[0]) - charge,weight=loadtxt(fileName+fileId+"/Results_extra",unpack=True,usecols=[0,2]) - - cond= (charge==0) - energy=energy[cond] - weight=weight[cond] - print "data shape:", shape(energy) - - Edata =log10(energy) - flux,ener=histogram(Edata,nbBins,range=(Emin,Emax),weights=weight) - ener = 10**ener - enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2 - binSize=ener[1:nbBins+1]-ener[0:nbBins] - flux = enercenter**2 * flux/binSize - - maxflux = shape(energy)[0] - - if argv[1] == "EGMF": - ax.plot(enercenter,flux[:]/maxflux,color=color[ind],drawstyle='steps-mid', - label="$J_\gamma$ $10^{-%.0f"%float(fileId)+"}$Gauss") - elif argv[1] == "EBL": - ax.plot(enercenter,flux[:]/maxflux,color=color[ind],drawstyle='steps-mid', - label=Model[int(fileId)-1]) - # Results from Elmag - Elmag=fileName+fileId+"/spec_diff_test" - if os.path.isfile(Elmag): - energy,fluxGamma,fluxLepton = loadtxt(Elmag,unpack=True,usecols=[0,1,2]) - energy=energy*10**(-9) - flux=fluxGamma+fluxLepton - flux=flux*2.5e-12 - ax.plot(energy,flux,"-"+color[ind],label="Elmag - "+Model[int(fileId)-1]) - - ind=ind+1 - -def fit(E_ic): - return me*1e3/4*sqrt(3e9/Ecmb)*1e-9*sqrt(E_ic) #GeV - -xfit = logspace(log10(10**(-3)),log10(2*10**4),200) -yfit = fit(xfit) -ax.plot(xfit,yfit,'--m',linewidth=2,label="$E^{0.5}$") - -ax.set_xscale('log') -ax.set_yscale('log') -ax.grid(b=True,which='major') -ax.legend(loc="best") -ax.set_xlabel("energy [GeV]") -ax.set_ylabel("$E^2 \\frac{dN}{dE}$ [GeV.photon]") - -plt.show() diff --git a/Time_delay.py b/Time_delay.py deleted file mode 100755 index 50ed2d2..0000000 --- a/Time_delay.py +++ /dev/null @@ -1,96 +0,0 @@ -#!/usr/bin/python - -from sys import argv -import os.path -import matplotlib.pyplot as plt -from Constantes import * -from scipy.integrate import quad - -if argv[1] == "EGMF": - fileName="Results_EGMF" -elif argv[1] == "EBL": - fileName="Results_EBL" - B=10**(-15) - Model=["Kneiske and Doll - 'best fit'","Kneiske and Doll - 'lower limit'", - "Franceschini","Finke and Al","Gilmore and Al","Dominguez and Al"] -elif argv[1] == "prec": - fileName="Results_prec" - B=10**(-17) -else: - print "Used: ./Energy_distribution.py arg1 ind1 ind2 ..." - print " arg1 = EGMF or EBL " - print " ind1 = file number \n" - print "Examples: " - print " to study EGMF14: ./Energy_distribution.py EGMF 14 " - print " to study EBL1 and EBL2: ./Energy_distribution.py EBL 1 2 \n" - exit() - -# figure: energy distribution -#============================= -zSource = 0.0308 -nbBins = 100 -convert= 180/pi -color=['b','r','g','c','m','y'] - -fig = plt.figure() -ax = fig.add_subplot(111) - -tlim = quad(comobileTime,zSource,0)[0]/yr -print tlim - -time= loadtxt(fileName+argv[2]+"/Results_position",unpack=True,usecols=[0]) -time=time[time>0] -time=time/tlim -tmax = max(time) -tmin = min(time) -print tmin, tmax - -ind = 0 -for fileId in argv[2:]: - time= loadtxt(fileName+fileId+"/Results_position",unpack=True,usecols=[0]) - charge,weight,gen=loadtxt(fileName+fileId+"/Results_extra",unpack=True,usecols=[0,2,3]) - time=time/tlim - - prop = float(shape(time[time<0])[0])/shape(time)[0] - print "data shape",fileName+fileId,":", shape(time), "negative time:", shape(time[time<0]), "~", prop - - cond= (gen<10) & (time>0) & (charge==0) - time=time[cond] - weight=weight[cond] - time=log10(time) - - dN,dt=histogram(time,nbBins,range=(log10(tmin),log10(tmax)),weights=weight) - #time=10**time - timecenter=(dt[1:nbBins+1]+dt[0:nbBins])/2 - binSize=dt[1:nbBins+1]-dt[0:nbBins] - dNdt=dN/binSize - - Nmax=shape(time)[0] - - if argv[1] == "EGMF": - #ax.hist(time,nbBins,weights=weight,range=(log10(tmin),log10(tmax)),log=1,facecolor=color[ind], - # alpha=.5, label="$10^{-%.0f"%float(fileId)+"}$Gauss") - ax.plot(timecenter[dNdt!=0],dNdt[dNdt!=0]/Nmax,color=color[ind],linestyle="steps-mid", - label="$10^{-%.0f"%float(fileId)+"}$Gauss") - elif argv[1] == "prec": - ax.hist(time,nbBins,weights=weight,range=(log10(tmin),log10(tmax)),log=1,facecolor=color[ind], - alpha=.5)#, label="prec = $10^{-%.0f"%float(fileId)+"}$") - ax.plot(timecenter[dNdt!=0],dNdt[dNdt!=0],"-"+color[ind], - label="prec = $10^{-%.0f"%float(fileId)+"}$") - elif argv[1] == "EBL": - #plt.hist(time,nbBins,weights=weight,range=(tmin,tmax),log=1,facecolor=color[ind],alpha=.5, - ax.plot(timecenter,dNdt,"."+color[ind], - label=Model[int(fileId)-1]) - - ind=ind+1 - -#ax.set_xscale('log') -ax.set_yscale('log') -ax.grid(b=True,which='major') -ax.legend(loc="best") -if argv[1]=="prec": - ax.legend(loc="best",title="$10^{-17}$Gauss") -ax.set_xlabel("Time [$\Delta t / t$ log scale]") -ax.set_ylabel("$t\ dN/dt$") - -plt.show() -- libgit2 0.21.2