From f42463907fb1bec9c6dc7d526d4c07d3220ead49 Mon Sep 17 00:00:00 2001 From: Thomas Fitoussi Date: Thu, 21 May 2015 10:44:11 +0200 Subject: [PATCH] Improve angular distribution (plot in log-log) --- Analysis.py | 2 +- Modules/Analytic.py | 9 ++++++++- Modules/Analytic.pyc | Bin 3203 -> 0 bytes Modules/Angle.py | 74 ++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------- Modules/Angle.pyc | Bin 4441 -> 0 bytes Modules/Constants.pyc | Bin 1107 -> 0 bytes Modules/Generation.py | 5 +++++ Modules/Generation.pyc | Bin 1137 -> 0 bytes Modules/Map.py | 29 ++++++++++++++++++++++------- Modules/Map.pyc | Bin 3505 -> 0 bytes Modules/Spectrum.py | 34 +++++++++++++++++++++++++++++----- Modules/Spectrum.pyc | Bin 4946 -> 0 bytes Modules/Timing.py | 150 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------- Modules/Timing.pyc | Bin 2166 -> 0 bytes Modules/Weight.py | 15 +++++++++++++++ Modules/Weight.pyc | Bin 0 -> 2434 bytes 16 files changed, 269 insertions(+), 49 deletions(-) create mode 100644 Modules/Weight.pyc diff --git a/Analysis.py b/Analysis.py index cf3eb73..22619ea 100755 --- a/Analysis.py +++ b/Analysis.py @@ -26,7 +26,7 @@ if shape(argv)[0] < 3: error("not enough argument") elif argv[1] == "spectrum": - drawSpectrum(argv[2:]) + drawSpectrum(argv[2:],PlotAnalytic=True) elif argv[1] == "spectrumLept": drawSpectrum_normaLepton(argv[2:]) diff --git a/Modules/Analytic.py b/Modules/Analytic.py index af5dfa4..03dbcc4 100644 --- a/Modules/Analytic.py +++ b/Modules/Analytic.py @@ -1,5 +1,5 @@ from Constants import * -from numpy import sqrt, abs, sin, arcsin +from numpy import sqrt, abs, cos, sin, arcsin from Read import ReadE_source, ReadD_source, ReadEGMF # Value for analytic expression @@ -21,6 +21,13 @@ def Analytic_theta(E_ic,fileId): delta=Dic0/(2*RL0)*((Esource/2)**2 /Ee2 -1) return abs(arcsin(delta_to_theta*sin(delta))*degre) +def Analytic_delay(theta, E0_source, distSource): + delta_ic = arcsin(lambda_gg(E0_source)/distSource*theta/degre) + #c_delta_t = lambda_gg(E0_source)*(1-cos(delta_ic))-distSource*(1-cos(theta)) + #return abs(c_delta_t*Mpc/c) + coef = lambda_gg(E0_source)/(2*c/Mpc) + return coef*delta_ic#**2 + # Compton accumulation def ECompton_threshold(Compton_threshold = 0.005): return Compton_threshold/(4/3*Ecmb/me*1e-3) *me*1e-6 #GeV diff --git a/Modules/Analytic.pyc b/Modules/Analytic.pyc index aebef08..9e0d698 100644 Binary files a/Modules/Analytic.pyc and b/Modules/Analytic.pyc differ diff --git a/Modules/Angle.py b/Modules/Angle.py index 35e305c..f5fe33c 100644 --- a/Modules/Angle.py +++ b/Modules/Angle.py @@ -3,12 +3,18 @@ from matplotlib.pyplot import figure, show, ylim, setp from matplotlib import gridspec from Read import ReadMomentum, ReadPosition, ReadEnergy, ReadWeight, ReadNphot_source from Analytic import Analytic_theta +from Map import thetaphi from Constants import degre Eliminf = 1e0 #GeV Elimsup = 1e5 #GeV def compute_theta(fileName): + ''' + Compute arrival angle + Input: directory name + Output: energy, weight, arrival angle + ''' position = ReadPosition(fileName) momentum = ReadMomentum(fileName) energy = ReadEnergy(fileName) @@ -29,12 +35,17 @@ def compute_theta(fileName): return energy, weight, theta def angle_vs_energy(fileName): + ''' + Compute arrival angle versus energy + Input: directory name + Output: energy, arrival angle + ''' energy,weight,theta = compute_theta(fileName) nbBins = 40 energy =log10(energy) - angle,ener =histogram(energy,nbBins,weights=theta) - nb,ener =histogram(energy,nbBins) + angle,ener =histogram(energy,nbBins,weights=weight*theta) + nb,ener =histogram(energy,nbBins,weights=weight) ener = 10**ener enercenter=(ener[1:nbBins+1]+ener[0:nbBins])/2 angle = angle/nb @@ -42,6 +53,11 @@ def angle_vs_energy(fileName): return enercenter, angle def drawAngle(files): + ''' + Plot arrival angle versus energy + analytic expression + Input: list of directories + Output: graph of arrival angle versus energy + ''' fig = figure() gs = gridspec.GridSpec(2, 1, height_ratios=[4,1]) fig.subplots_adjust(hspace=0.001) @@ -75,31 +91,49 @@ def drawAngle(files): 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 - maxflux = ReadNphot_source(fileName) - dndtheta2 /= maxflux - - return thetacenter, dndtheta2 - -def drawRadial(files): +############################################################################################## +def radial(fileId,thetamax=90): + ''' + Compute flux versus arrival angle + Input: directory name + Input (optionnal): maximal angle desired (by default full sky) + Output: arrival angle, flux (dn/dtheta) + ''' + nbBins = 50 + energy,weight,theta = compute_theta(fileId) + + theta=log10(theta) + dn,angle = histogram(theta,nbBins,range=[log10(1e-6),log10(thetamax)],weights=weight) + angle=10**angle + thetacenter = (angle[1:nbBins+1]+angle[0:nbBins])/2 + dtheta = angle[1:nbBins+1]-angle[0:nbBins] + dndtheta = dn/dtheta + maxflux = ReadNphot_source(fileId) + dndtheta /= maxflux + + print "integral: ", sum(dtheta*dndtheta) + + return thetacenter, dndtheta + +def drawRadial(files,thetamax=90): + ''' + Plot flux versus arrival angle + Input: list of directories + Input (optionnal): maximal angle desired (by default full sky) + Output: graph of flux versus angle + ''' fig = figure() ax = fig.add_subplot(111) for fileId in files: - theta,dndtheta2=radial(fileId) + theta,dndtheta2=radial(fileId,thetamax) ax.plot(theta,dndtheta2,drawstyle='steps-mid',label=fileId) - ax.legend(loc="best",title="%3d - %3d GeV"%(Eliminf,Elimsup)) + ax.legend(loc="best") + 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$ [deg$^{-2}$]") + ax.set_xlabel("$\\theta$ [deg]") + ax.set_ylabel("$dN/d\\theta$ [deg$^{-1}$]") show() diff --git a/Modules/Angle.pyc b/Modules/Angle.pyc index 694484e..8e9cb79 100644 Binary files a/Modules/Angle.pyc and b/Modules/Angle.pyc differ diff --git a/Modules/Constants.pyc b/Modules/Constants.pyc index 88ddc03..fe7442f 100644 Binary files a/Modules/Constants.pyc and b/Modules/Constants.pyc differ diff --git a/Modules/Generation.py b/Modules/Generation.py index b038938..cea5271 100644 --- a/Modules/Generation.py +++ b/Modules/Generation.py @@ -4,6 +4,11 @@ from Read import ReadGeneration, ReadWeight def drawGeneration(files): + ''' + Plot histogram of particles'generations + Input: list of directories + Output: histogram of generations + ''' fig = figure() ax = fig.add_subplot(111) nbBins = 7 diff --git a/Modules/Generation.pyc b/Modules/Generation.pyc index 79de4ef..257b80e 100644 Binary files a/Modules/Generation.pyc and b/Modules/Generation.pyc differ diff --git a/Modules/Map.py b/Modules/Map.py index ddeeb2c..332a8f5 100644 --- a/Modules/Map.py +++ b/Modules/Map.py @@ -1,4 +1,4 @@ -from numpy import shape, pi, histogram2d, arange, newaxis, exp, log, array +from numpy import shape, pi, histogram2d, arange, newaxis, exp, log, array, sin, cos from scipy.signal import convolve2d from matplotlib.pyplot import figure, show, matshow from matplotlib.colors import LogNorm @@ -7,14 +7,21 @@ from Read import ReadD_source, Readz_source, ReadEGMF from Constants import degre thetalim = 180 -philim = 90 -step = 1 #deg +philim = 180 +step = 0.25 #deg nbBins = array([int(2*philim/step)+1,int(2*thetalim/step)+1]) philim=float(philim) thetalim=float(thetalim) limits = [[-thetalim,thetalim],[-philim,philim]] def thetaphi(fileId,Elim=0.1): + ''' + Compute arrival angles (theta,phi) selected with energy upper than Elim + normalize to 1 source photon + Input: directory name + Input (optional): lower energy limit in GeV (default = 100MeV) + Output: energy (GeV), theta, phi (degre), weight + ''' energy = ReadEnergy(fileId) weight = ReadWeight(fileId) thetaDir,phiDir = ReadMomentumAngle(fileId)*degre @@ -42,7 +49,15 @@ def Gaussian2D(center=[0,0],fwhm=0.01): y = y[:,newaxis] return exp(-4*log(2) * ((x-center[0])**2 + (y+center[1])**2) / fwhm**2) -def drawMap(files,Elim=0.1,source=isotrop(),borne=[180,90]): +def drawMap(files,Elim=0.1,source=isotrop(),borne=[180,180]): + ''' + Plot 2D histogram of arrival direction of the photons with energy upper than Elim + Input: list of directories + Input (optional): lower energy limit (in GeV (default = 100MeV) + Input (optional): nature of the source (default isotrop) + Input (optional): limits on theta, phi + Output: 2D histogram of theta cos(phi), theta sin(phi) + ''' print "Cut at", Elim, "GeV" fig = figure() nbplots = len(files) @@ -51,7 +66,7 @@ def drawMap(files,Elim=0.1,source=isotrop(),borne=[180,90]): for fileId in files: energy, theta, phi, weight = thetaphi(fileId,Elim) - H,xedges,yedges = histogram2d(phi,theta,nbBins,weights=weight,range=limits) + H,xedges,yedges = histogram2d(theta*cos(phi),theta*sin(phi),nbBins,weights=weight,range=limits) conv=convolve2d(H,source,boundary="wrap",mode="same") ax = fig.add_subplot(100+nbplots*10+ind) @@ -62,8 +77,8 @@ def drawMap(files,Elim=0.1,source=isotrop(),borne=[180,90]): cbar.ax.set_ylabel("counts normalize to 1 photon emitted") ax.set_xlim((-borne[0],borne[0])) ax.set_ylim((-borne[1],borne[1])) - ax.set_xlabel("$\\theta$ [deg]") - ax.set_ylabel("$\\phi$ [deg]") + ax.set_xlabel("$\\theta$ $\\cos(\\phi)$ [deg]") + ax.set_ylabel("$\\theta$ $\\sin(\\phi)$ [deg]") ax.grid(b=True,which='major') zsource=Readz_source(fileId) Dsource=ReadD_source(fileId) diff --git a/Modules/Map.pyc b/Modules/Map.pyc index 42389d4..6d359bc 100644 Binary files a/Modules/Map.pyc and b/Modules/Map.pyc differ diff --git a/Modules/Spectrum.py b/Modules/Spectrum.py index 5eb7603..76ad36b 100644 --- a/Modules/Spectrum.py +++ b/Modules/Spectrum.py @@ -7,6 +7,11 @@ from Read import ReadNphot_source, ReadGeneration from Analytic import Analytic_flux, Ethreshold_gg def spectrum(fileId,gen=-1): + ''' + Compute flux versus energy + Input: directory name, generation (optional, by default: all generation) + Output: energy, flux + ''' energy = ReadEnergy(fileId) weight = ReadWeight(fileId) @@ -30,12 +35,19 @@ def spectrum(fileId,gen=-1): return enercenter,flux ########################################################################################### -def drawSpectrum(files,PlotElmag=False,PlotAnalytic=False): +def drawSpectrum(files,gen=-1,PlotElmag=False,PlotAnalytic=False): + ''' + Plot flux versus energy + Input: list of directories + Input (optional, default=all): selection on the generation + Input (optional, default=False): Plot Elmag spectrum if exists, plot analytic expression + Output: graph spectrum normalize to 1 source photon + ''' fig = figure() ax = fig.add_subplot(111) for fileId in files: - energy,flux = spectrum(fileId) + energy,flux = spectrum(fileId,gen) if PlotElmag: maxflux = max(flux)*(.8) p=ax.plot(energy,flux,drawstyle='steps-mid',label=fileId) @@ -47,9 +59,10 @@ def drawSpectrum(files,PlotElmag=False,PlotAnalytic=False): ax.plot(energy,flux,"--",color=p[0].get_color(),label="Elmag - "+fileId) if PlotAnalytic: - alpha=2.8e3 + minEner=min(energy) + alpha=flux[energy==minEner]/sqrt(minEner) ax.plot(energy,alpha*sqrt(energy),'--m',linewidth=2) - ax.plot(energy,energy**(0)*max(flux),'--m',linewidth=2) + ax.plot(energy,energy**(0)*max(flux)*0.95,'--m',linewidth=2) ax.axvline(x=Ethreshold_gg(), ymin=0., ymax = 1., color='r', linewidth=2) ax.set_xscale('log') @@ -65,6 +78,11 @@ def drawSpectrum(files,PlotElmag=False,PlotAnalytic=False): show() def drawSpectrum_normaLepton(files): + ''' + Plot flux versus energy compare with analytic expression + Input: list of directories + Output: graph spectrum normalize to 1 lepton generated + ''' fig = figure() gs = gridspec.GridSpec(2, 1, height_ratios=[4,1]) fig.subplots_adjust(hspace=0.001) @@ -101,6 +119,11 @@ def drawSpectrum_normaLepton(files): ######################################################################################## def drawSpectrumGen(fileId): + ''' + Plot flux versus energy generation by generation + Input: directory name + Output: graph spectrum normalize to 1 photon source + ''' fig = figure() ax = fig.add_subplot(111) @@ -110,8 +133,9 @@ def drawSpectrumGen(fileId): energy,flux = spectrum(fileId) p=ax.plot(energy,flux,drawstyle='steps-mid',color='k',label="gen = all") - alpha=2e3 + minEner=min(energy) + alpha=flux[energy==minEner]/sqrt(minEner) ax.plot(energy,sqrt(energy)*alpha,'--m',linewidth=2) ax.plot(energy,energy**(0)*max(flux),'--m',linewidth=2) ax.axvline(x=Ethreshold_gg(), ymin=0., ymax = 1., color='r', linewidth=2) diff --git a/Modules/Spectrum.pyc b/Modules/Spectrum.pyc index 69b8fd5..a95375e 100644 Binary files a/Modules/Spectrum.pyc and b/Modules/Spectrum.pyc differ diff --git a/Modules/Timing.py b/Modules/Timing.py index b994ed0..472205d 100644 --- a/Modules/Timing.py +++ b/Modules/Timing.py @@ -1,45 +1,165 @@ -from numpy import histogram, log10, shape +from numpy import histogram, log10, shape, sqrt, arccos from matplotlib.pyplot import figure, show from scipy.integrate import quad -from Read import ReadTime, ReadWeight, Readz_source, ReadNphot_source -from Constants import yr +from Read import ReadTime, ReadWeight, ReadGeneration, ReadE_source, ReadD_source +from Read import ReadNphot_source,ReadPosition,ReadMomentum,ReadNbLeptonsProd +from Constants import yr, degre from Integrand import comobileTime +from Analytic import Analytic_delay -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])[0], "~", prop +def timing(fileId,gen=-1): + ''' + Compute flux versus time delay + Input: directory name + Input (optional): generation (default = all) + Output: time delay (sec), flux + ''' + time=ReadTime(fileId) + weight=ReadWeight(fileId) + + if gen != -1: + generation = ReadGeneration(fileId) + time=time[generation==gen] + weight=weight[generation==gen] + prop = float(shape(time[time<0])[0])/shape(time)[0] + print "gen", int(gen), "->", shape(time)[0], "events", "negative time:",shape(time[time<0])[0], "~", prop + else: + prop = float(shape(time[time<0])[0])/shape(time)[0] + print "file", fileId, "->", shape(time)[0], "events", "negative time:",shape(time[time<0])[0], "~", prop time=time[time>0] weight=weight[time>0] nbBins = 100 - zSource = Readz_source(fileName) - tlim = quad(comobileTime,zSource,0)[0]/yr - time = log10(time/tlim) + time = log10(time*yr) 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=ReadNphot_source(fileName) # Nb photons emitted by the source + Nmax=ReadNphot_source(fileId) # Nb photons emitted by the source dNdt=dN/(Nmax*binSize) - maxflux = ReadNphot_source(fileName) + maxflux = ReadNphot_source(fileId) dNdt /= maxflux - return timecenter, dNdt + return 10**timecenter, dNdt def drawTiming(files): + ''' + Plot flux versus time delay + Input: list of directories + Output: graph of flux versus time delay + ''' 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_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("$t\ dN/dt$") + + show() + +def drawTimingGen(fileId): + fig = figure() + ax = fig.add_subplot(111) + for gen in list(set(ReadGeneration(fileId))): + time,dNdt = timing(fileId,gen) + ax.plot(time[dNdt!=0],dNdt[dNdt!=0],linestyle="steps-mid",label="gen="+str(int(gen))) + + time,dNdt = timing(fileId) + ax.plot(time[dNdt!=0],dNdt[dNdt!=0],color='k',linestyle="steps-mid",label="gen=all") + + ax.set_xscale('log') 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_xlabel("Time delay [s]") ax.set_ylabel("$t\ dN/dt$") show() + +######## delay versus angle ############### +def delay_vs_angle(fileId,gen=-1): + ''' + Compute time delay versus arrival angle + Input: directory name + Input (optional): generation (default = all) + Output: arrival angle [degre], time delay [sec] + ''' + position = ReadPosition(fileId) + momentum = ReadMomentum(fileId) + weight = ReadWeight(fileId) + time=ReadTime(fileId) + + 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] + + if gen != -1 : + generation = ReadGeneration(fileId) + cond=(abs(costheta)<=1) & (generation==gen) & (time>0) + else: + cond=(abs(costheta)<=1) & (time>0) + + prop = float(shape(time[time<0])[0])/shape(time)[0] + theta = arccos(costheta[cond])*degre + weight= weight[cond] + time=time[cond]*yr + print "gen", int(gen), "->", shape(time)[0], "events", "negative time:",shape(time[time<0])[0], "~", prop + cond=(theta!=0) + return theta[cond], time[cond], weight[cond] + +def delay_vs_angle_histo(fileId,gen=-1): + theta, time, weight = delay_vs_angle(fileId,gen) + nbBins = 100 + + theta=log10(theta) + dt,dtheta=histogram(theta,nbBins,weights=time) + dN,dtheta=histogram(theta,nbBins,weights=weight) + dtheta=10**dtheta + thetacenter=(dtheta[1:nbBins+1]+dtheta[0:nbBins])/2 + timing = dt/dN + + return thetacenter, timing + +def drawDelay_vs_angle(fileId): + ''' + Plot angle versus time delay, generation by generation + Input: directory name + Output: graph of angle versus time delay, generation by generation + ''' + fig = figure() + ax = fig.add_subplot(111) + + Nmax = ReadNphot_source(fileId) + + for gen in list(set(ReadGeneration(fileId))): + angle,delay,weight = delay_vs_angle(fileId,gen) + ax.plot(angle,delay/weight/Nmax,".",label="gen="+str(int(gen))) + + #angle,delay = delay_vs_angle_histo(fileId,gen) + #ax.plot(angle,delay,linestyle="steps-mid",label="gen="+str(int(gen))) + + angle,delay = delay_vs_angle_histo(fileId) + ax.plot(angle,delay/Nmax,color='k',linestyle="steps-mid",label="gen=all") + + Esource = ReadE_source(fileId) + dSource = ReadD_source(fileId) + yfit = Analytic_delay(angle,Esource,dSource) + ax.plot(angle,yfit,'--k',linewidth=2) + + ax.set_xscale('log') + ax.set_yscale('log') + ax.grid(b=True,which='major') + ax.legend(loc="best") + ax.set_xlabel("$\\theta $[deg]") + ax.set_ylabel("Time delay [s]") + + show() diff --git a/Modules/Timing.pyc b/Modules/Timing.pyc index 8562513..50eeae8 100644 Binary files a/Modules/Timing.pyc and b/Modules/Timing.pyc differ diff --git a/Modules/Weight.py b/Modules/Weight.py index 1a94bf1..d0e9a6e 100644 --- a/Modules/Weight.py +++ b/Modules/Weight.py @@ -4,6 +4,11 @@ from Read import ReadGeneration, ReadWeight, ReadEnergy def drawWeightHisto(fileId): + ''' + Plot particles'weights histogram, generation by generation + Input: directory name + Output: particles'weights histogram + ''' fig = figure() ax = fig.add_subplot(111) nbBins =100 @@ -24,6 +29,11 @@ def drawWeightHisto(fileId): ####################################################################################### def weightVsEnergy(fileId,gen=-1): + ''' + Compute particles'weight versus energy + Input: directory name + Input (optional): generation (default = all) + ''' energy = ReadEnergy(fileId) weight = ReadWeight(fileId) nbBins = 100 @@ -43,6 +53,11 @@ def weightVsEnergy(fileId,gen=-1): return enercenter, weight def drawWeightVsEnergy(fileId): + ''' + Plot particles'weight versus energy, generation by generation + Input: directory name + Output: graph of particles'weight versus energy + ''' fig = figure() ax = fig.add_subplot(111) diff --git a/Modules/Weight.pyc b/Modules/Weight.pyc new file mode 100644 index 0000000..576c2b2 Binary files /dev/null and b/Modules/Weight.pyc differ -- libgit2 0.21.2