import matplotlib.cm as colormap from matplotlib.pyplot import figure, show, legend, gca, setp from matplotlib.gridspec import GridSpec from numpy import sort, log10, histogram, pi, arange, where, ma, zeros, append from numpy import rot90, flipud, histogram2d, logspace, linspace, size, loadtxt from numpy import array, argsort, nditer, mean from mpl_toolkits.mplot3d import Axes3D from src.read import select_events, ReadProfile from src.analytic import degre, yr, Analytic_delay_vs_theta, Analytic_observables_vs_energy import matplotlib as mpl label_size = 20 mpl.rcParams['xtick.labelsize'] = label_size mpl.rcParams['ytick.labelsize'] = label_size def weighted_median(values, weights): ''' compute the weighted median of values list. The weighted median is computed as follows: 1- sort both lists (values and weights) based on values. 2- select the 0.5 point from the weights and return the corresponding values as results e.g. values = [1, 3, 0] and weights=[0.1, 0.3, 0.6] assuming weights are probabilities. sorted values = [0, 1, 3] and corresponding sorted weights = [0.6, 0.1, 0.3] the 0.5 point on weight corresponds to the first item which is 0. so the weighted median is 0.''' #convert the weights into probabilities sum_weights = sum(weights) weights = array([(w*1.0)/sum_weights for w in weights]) #sort values and weights based on values values = array(values) sorted_indices = argsort(values) values_sorted = values[sorted_indices] weights_sorted = weights[sorted_indices] #select the median point it = nditer(weights_sorted, flags=['f_index']) accumulative_probability = 0 median_index = -1 while not it.finished: accumulative_probability += it[0] if accumulative_probability > 0.5: median_index = it.index return values_sorted[median_index] elif accumulative_probability == 0.5: median_index = it.index it.iternext() next_median_index = it.index return mean(values_sorted[[median_index, next_median_index]]) it.iternext() return values_sorted[median_index] #=========================== CORELATION BETWEEN OBSERVABLES ===========================# def observables_vs_energy(energy,theta,delay,weight,Erange=[1e-3,1e5],median=False,nbBins=-1):#28): # nbBins=16, Erange=[1e-1,1e3] <=> Arlen 2014 (energy binning) if (nbBins==-1): bins=append(logspace(log10(Erange[0]),log10(0.9),100),logspace(0,log10(Erange[1]),16)) nbBins=size(bins)-1 if (nbBins==-2): # <=> Arlen 2014 (energy binning) nbBins=16 Erange=[1e-1,1e3] bins=logspace(log10(Erange[0]),log10(Erange[1]),nbBins+1) else: bins=logspace(log10(Erange[0]),log10(Erange[1]),nbBins+1) ener = (bins[1:nbBins+1]*bins[0:nbBins])**0.5 #dtheta, dE = histogram(energy,bins,weights=weight*theta) #dt, dE = histogram(energy,bins,weights=weight*delay) #dN, dE = histogram(energy,bins,weights=weight) dt = zeros(nbBins,dtype="float64") dtheta = zeros(nbBins,dtype="float64") if median==False: for i in arange(nbBins): mask = (energy >= bins[i]) & (energy < bins[i+1]) w = sum(weight[mask]) t = sum(weight[mask]*delay[mask]) angle = sum(weight[mask]*theta[mask]) if w !=0: dt[i]=t/w dtheta[i]=angle/w else: for i in arange(nbBins): mask = (energy >= bins[i]) & (energy < bins[i+1]) if True in mask: dt[i] = weighted_median(delay[mask],weight[mask]) dtheta[i] = weighted_median(theta[mask],weight[mask]) else: dt[i] = 0 dtheta[i] = 0 return ener, dtheta, dt def observables_vs_delay(energy,theta,delay,weight,time_range=[10**-0.5,10**8.5],median=False,nbBins=100):#9): # <=> Taylor 2011 (time binning) bins=logspace(log10(time_range[0]),log10(time_range[1]),nbBins+1) time = (bins[1:nbBins+1]*bins[0:nbBins])**0.5 dE = zeros(nbBins,dtype="float64") dtheta = zeros(nbBins,dtype="float64") if median==False: for i in arange(nbBins): mask = (delay >= bins[i]) & (delay < bins[i+1]) w = sum(weight[mask]) E = sum(weight[mask]*energy[mask]) angle = sum(weight[mask]*theta[mask]) if w !=0: dE[i]=E/w dtheta[i]=angle/w else: for i in arange(nbBins): mask = (delay >= bins[i]) & (delay < bins[i+1]) if True in mask: dE[i] = weighted_median(energy[mask],weight[mask]) dtheta[i] = weighted_median(theta[mask],weight[mask]) else: dE[i] = 0 dtheta[i] = 0 return time, dtheta, dE def drawObservables(fileId,psf=180,Nb=28,plot_generation_density=False,plot_others_codes=False,one_figure=True): #ax = figure(figsize=(12,9)).add_subplot(111,projection='3d') if one_figure: fig = figure(figsize=(20,18)) gs = GridSpec(2, 2, height_ratios=[1,1], width_ratios=[1,1]) fig.subplots_adjust(hspace=0,wspace=0) ax0 = fig.add_subplot(gs[0]) ax2 = fig.add_subplot(gs[2],sharex=ax0) ax3 = fig.add_subplot(gs[3],sharey=ax2) else: fig0 = figure(figsize=(12,9)) ax0 = fig0.add_subplot(111) fig2 = figure(figsize=(12,9)) ax2 = fig2.add_subplot(111) fig3 = figure(figsize=(12,9)) ax3 = fig3.add_subplot(111) fileId = "Simulations/"+fileId i=0 nbBins = 250 theta_range = [1e-5,psf] delay_range = [1e-4,5e9] energy_range = [1e-3,1e4] weight,energy,delay,theta,phi,Esource,generation = select_events(fileId,Erange=energy_range) delay /= yr theta /= degre print "theta range: [",min(theta),", ",max(theta),"] (degre)" print "delay range: [",min(delay),", ",max(delay),"] (yr)" print "energy range: [",min(energy),", ",max(energy),"] (GeV)" print "weight range: [",min(weight),", ",max(weight),"]" colors=['b','g'] cond = (theta>0) & (delay>0) if plot_generation_density: cmaps=[colormap.Blues,colormap.Greens,colormap.Reds] for gen in [2,4]:#sort(list(set(generation))): cond = cond & (generation==gen) H, xedges, yedges = histogram2d(log10(energy[cond]),log10(delay[cond]),bins=nbBins,weights=weight[cond]) H = flipud(rot90(ma.masked_where(H==0,H))) im1 = ax0.pcolormesh(10**xedges,10**yedges,log10(H),cmap=cmaps[i]) H, xedges, yedges = histogram2d(log10(energy[cond]),log10(theta[cond]),bins=nbBins,weights=weight[cond]) H = flipud(rot90(ma.masked_where(H==0,H))) im2 = ax2.pcolormesh(10**xedges,10**yedges,log10(H),cmap=cmaps[i]) H, xedges, yedges = histogram2d(log10(delay[cond]),log10(theta[cond]),bins=nbBins,weights=weight[cond]) H = flipud(rot90(ma.masked_where(H==0,H))) im3 = ax3.pcolormesh(10**xedges,10**yedges,log10(H),cmap=cmaps[i]) ener,angle,dt = observables_vs_energy(energy[cond],theta[cond],delay[cond],weight[cond],nbBins=Nb,median=False) ax0.plot(ener,dt,color=colors[i],linewidth=2,label="gen = %.0f"%(i+1)) ax2.plot(ener,angle,color=colors[i],linewidth=2,label="gen = %.0f"%(i+1)) dt,angle,ener = observables_vs_delay(energy[cond],theta[cond],delay[cond],weight[cond],nbBins=Nb,median=False) ax3.plot(dt,angle,color=colors[i],linewidth=2,label="gen = %.0f"%(i+1)) i+=1 else: H, xedges, yedges = histogram2d(log10(energy[cond]),log10(delay[cond]),bins=nbBins,weights=weight[cond]) H = flipud(rot90(ma.masked_where(H==0,H))) im1 = ax0.pcolormesh(10**xedges,10**yedges,log10(H),cmap=colormap.YlOrBr) H, xedges, yedges = histogram2d(log10(energy[cond]),log10(theta[cond]),bins=nbBins,weights=weight[cond]) H = flipud(rot90(ma.masked_where(H==0,H))) im2 = ax2.pcolormesh(10**xedges,10**yedges,log10(H),cmap=colormap.YlOrBr) H, xedges, yedges = histogram2d(log10(delay[cond]),log10(theta[cond]),bins=nbBins,weights=weight[cond]) H = flipud(rot90(ma.masked_where(H==0,H))) im3 = ax3.pcolormesh(10**xedges,10**yedges,log10(H),cmap=colormap.YlOrBr) for gen in [2,4]:#sort(list(set(generation))): cond = (generation==gen) ener,angle,dt = observables_vs_energy(energy[cond],theta[cond],delay[cond],weight[cond],nbBins=Nb,median=False) ax0.plot(ener,dt,color=colors[i],linewidth=2,label="gen = %.0f"%(i+1)) ax2.plot(ener,angle,color=colors[i],linewidth=2,label="gen = %.0f"%(i+1)) dt,angle,ener = observables_vs_delay(energy[cond],theta[cond],delay[cond],weight[cond],nbBins=Nb,median=False) ax3.plot(dt,angle,color=colors[i],linewidth=2,label="gen = %.0f"%(i+1)) i+=1 ener,angle,dt = observables_vs_energy(energy,theta,delay,weight,nbBins=Nb,median=False) ax0.plot(ener,dt,color='k',linewidth=2,label="all gen") ax2.plot(ener,angle,color='k',linewidth=2,label="all gen") dt,angle,ener = observables_vs_delay(energy,theta,delay,weight,nbBins=Nb,median=False) ax3.plot(dt,angle,color='k',linewidth=2,label="all gen") nE = 5000 Emin = 1.e-3 Emax = 1e4 E = Emin*(Emax/Emin)**(arange(nE)/(nE-1.)) theta_fit, delay_fit, theta_fit2, delay_fit2 = Analytic_observables_vs_energy(E,"../"+fileId) ax0.plot(E,delay_fit/yr,'--b',linewidth=2) ax0.plot(E,delay_fit2/yr,'--g',linewidth=2) ax2.plot(E,theta_fit,'--b',linewidth=2) ax2.plot(E,theta_fit2,'--g',linewidth=2) ax3.plot(delay_fit[delay_fit.argsort()]/yr,theta_fit[delay_fit.argsort()],'--b',linewidth=2) ax3.plot(delay_fit2[delay_fit2.argsort()]/yr,theta_fit2[delay_fit2.argsort()],'--g',linewidth=2) #ax3.scatter(delay_fit/yr,theta_fit,color='b',marker="+") if plot_others_codes: # Results from Arlen 2014 - fig. 2a, 2b # ===================================== # - binning in energy : log, 16 bins between 1e-1 et 1e3 GeV # - PSF = 10 deg data = loadtxt(fileId+'/Arlen2014-fig2a.csv', delimiter=',') ax0.plot(data[:,0],data[:,1],color="m",linestyle='--',linewidth=2,label="Arlen 2014 - no EBL") data = loadtxt(fileId+'/Arlen2014-fig2b.csv', delimiter=',') ax0.plot(data[:,0],data[:,1],color="m",linestyle='-.',linewidth=2,label="Arlen 2014 - EBL") # Results from Taylor 2011 - fig. 2 # ================================= data = loadtxt(fileId+'/Taylor2011-fig2.csv', delimiter=',') ax0.plot(data[:,0]*1e-9,data[:,1],color="m",linestyle=':',linewidth=2,label="Taylor 2011") #ax.legend(loc="best") #ax.set_xlabel("$E$ [GeV]") #ax.set_ylabel("$\\theta$ [deg]") #ax.set_zlabel("$\\Delta t$ [s]") #ax3.set_title(fileId+" - selection %0.1f"%psf) if one_figure: ax0.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.) ax0.grid(b=True,which='major') ax0.set_xscale('log') ax0.set_yscale('log') ax0.set_xlim(energy_range) ax0.set_ylim(delay_range) ax0.set_ylabel("$\\Delta t$ [yrs]",fontsize=label_size) if one_figure: setp(ax0.get_xticklabels(), visible=False) else: ax0.set_xlabel("$E$ [GeV]",fontsize=label_size) ax0.legend(loc="best",fontsize=label_size) cbar1=fig0.colorbar(im1, ax=ax0) cbar1.ax.set_ylabel("counts [log]",fontsize=label_size) ax2.grid(b=True,which='major') ax2.set_xscale('log') ax2.set_yscale('log') ax2.set_xlim(energy_range) ax2.set_ylim(theta_range) ax2.set_xlabel("$E$ [GeV]",fontsize=label_size) ax2.set_ylabel("$\\theta$ [deg]",fontsize=label_size) if not one_figure: ax2.legend(loc="best",fontsize=label_size) cbar2=fig2.colorbar(im2, ax=ax2) cbar2.ax.set_ylabel("counts [log]",fontsize=label_size) ax3.grid(b=True,which='major') ax3.set_xscale('log') ax3.set_yscale('log') ax3.set_xlim(delay_range) ax3.set_ylim(theta_range) ax3.set_xlabel("$\\Delta t$ [yrs]",fontsize=label_size) if one_figure: setp(ax3.get_yticklabels(), visible=False) else: ax3.set_ylabel("$\\theta$ [deg]",fontsize=label_size) ax3.legend(loc="best",fontsize=label_size) cbar3=fig3.colorbar(im3, ax=ax3) cbar3.ax.set_ylabel("counts [log]",fontsize=label_size) show() def drawDelays_vs_energy(files,psf=180,plot_others_codes=False): ax1 = figure(figsize=(12,9)).add_subplot(111) nbBins = 100 for fileId0 in files: fileId = "Simulations/simple case/"+fileId0 weight, energy, delay, arrival_angle , theta, phi, gen = select_events(fileId) #dt,dtheta,ener = observables_vs_delay(energy,arrival_angle/degre,delay/yr,weight) #ax1.plot(ener,dt,marker='*',linestyle=":",linewidth=2) cond = (arrival_angle/degre < psf) & (gen%2==0) #PSF_Taylor2011(energy) energy = energy[cond] arrival_angle = arrival_angle[cond] delay = delay[cond] weight = weight[cond] ener,dtheta,dt = observables_vs_energy(energy,arrival_angle/degre,delay/yr,weight,nbBins=-2) p0, =ax1.plot(ener,dt,marker='*',linewidth=2,label=fileId0) nth = 5000 thmin = 0.1 thmax = 1e4 E = thmin*(thmax/thmin)**(arange(nth)/(nth-1.)) #ax1.plot(E,Analytic_observables_vs_energy(E,fileId)[1]/yr,color=p0.get_color(),linestyle="-") if plot_others_codes: # Results from Arlen 2014 - fig. 2a, 2b # ===================================== # - binning in energy : log, 16 bins between 1e-1 et 1e3 GeV # - PSF = 10 deg data = loadtxt(fileId+'/Arlen2014-fig2a.csv', delimiter=',') p1,= ax1.plot(data[:,0],data[:,1],color="r",linestyle='--',linewidth=2,label="Arlen 2014 - fig 2a") data = loadtxt(fileId+'/Arlen2014-fig2b.csv', delimiter=',') p2,= ax1.plot(data[:,0],data[:,1],color="r",linestyle='-.',linewidth=2,label="Arlen 2014 - fig 2b") # Results from Taylor 2011 - fig. 2 # ================================= #data = loadtxt(fileId+'/Taylor2011-fig2.csv', delimiter=',') #p3,= ax1.plot(data[:,0]*1e-9,data[:,1],color=p0.get_color(),linestyle=':',linewidth=2) #leg1 = legend([p0,p1,p2], ["our code","Arlen 2014 - fig 2a", # "Arlen 2014 - fig 2b"], loc=1) ax1.grid(b=True,which='major') ax1.set_xscale('log') ax1.set_yscale('log') #ax1.set_ylim([1,1e8]) #ax1.set_xlim([1e-1,1e3]) ax1.set_xlabel("$E$ [GeV]",fontsize=label_size) ax1.set_ylabel("$\\Delta t$ [yrs]",fontsize=label_size) ax1.legend(loc="best",fontsize=label_size)#2) #ax1 = gca().add_artist(leg1) show()