added updated cross-correlation function

This commit is contained in:
P.M. Sutter 2014-05-12 18:26:03 -04:00
parent fc19453800
commit 45a2b7f82f
2 changed files with 288 additions and 367 deletions

View file

@ -1,187 +1,58 @@
#+
# VIDE -- Void IDentification and Examination -- ./analysis/xcor.py
# Copyright (C) 2010-2014 Guilhem Lavaux
# Copyright (C) 2011-2014 P. M. Sutter
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; version 2 of the License.
#
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#+
#!/usr/bin/env python
#+
# VIDE -- Void IDentification and Examination -- ./crossCompare/analysis/mergerTree.py
# Copyright (C) 2010-2013 Guilhem Lavaux
# Copyright (C) 2011-2013 P. M. Sutter
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; version 2 of the License.
#
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#+
from void_python_tools.backend import *
import imp
import pickle
import argparse
import os
import string
import numpy as np import numpy as np
import matplotlib as mpl import matplotlib as mpl
mpl.use('Agg') mpl.use('Agg')
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import matplotlib.cm as cm import matplotlib.cm as cm
from matplotlib import rc from matplotlib import rc
from matplotlib.ticker import NullFormatter import xcorlib
import random
import sys
__all__=['computeCrossCor',]
# ------------------------------------------------------------------------------ def computeXcor(catalog,
figDir="./",
Nmesh = 256,
Nbin = 100
):
def computeCrossCor(catalogDir, # Computes and plots void-void and void-matter(galaxy) correlations
outputDir="./", logDir="./", # catalog: catalog to analyze
matchMethod="useID", dataPortion="central", # figDir: where to place plots
strictMatch=True, # Nmesh: number of grid cells in power spectrum calculation
pathToCTools="../../../c_tools"): # Nbin: number of grid cells in plot
# Computes void-void and void-matter(galaxy) correlations # Parameters
# baseCatalogDir: directory of catalog Lbox = catalog.boxLen # Boxlength
# compareCatagDir: directory of catalog 2 Lboxcut = 0.
# outputDir: directory to place outputs Lbox -= 2*Lboxcut
# logDir: directory to place log files
# matchMethod: "useID" to use unique IDs, "prox" to use overlap of Voronoi cells
# dataPortion: "central" or "all"
# strictMatch: if True, only attempt to match to trimmed catalog
# pathToCTools: path to location of VIDE c_tools directory
if not os.access(outputDir, os.F_OK): # Input particle arrays of shape (N,3)
os.makedirs(outputDir) xm = catalog.partPos # Halos / Galaxies / Dark matter
xv = catalog.voids[:].barycenter # Voids
with open(catalogDir+"/sample_info.dat", 'rb') as input:
sample = pickle.load(input)
print " Working with", sample.fullName, "...",
sys.stdout.flush()
sampleName = sample.fullName
# Sim parameters
Lbox = sample.boxLen # Boxlength [h^(-1)Mpc]
Lbox -= 2*Lboxcut # Reduced boxlength [h^(-1)Mpc]
Om = sample.omegaM # Omega_m
Ol = 1.-Om # Omega_l
z = sample.zRange[0] # Redshift
a = 1./(1.+z) # Scale factor
rho_m = Mpart*(Ni/Lbox)**3 # Background density [(h/Mpc)^3]
# Input files
voidDir = voidBaseDir+'/'+sampleDir+'/'
voidFilename1 = 'centers_central_'+sample.fullName+'.out'
voidFilename2 = 'voidDesc_central_'+sample.fullName+'.out'
# Read files
matter_file = open(matterDir+matterFilename,'r')
matter_data = []
count = 0
sep = 1e8
for (i,line) in enumerate(matter_file):
if i/int(sep) == i/sep: print str(round(i*100./Ni**3,1))+' percent of all particles read'
if random.random() > ss: continue
count += 1
matter_data.append(line.split(',')[0:3])
print str(count)+' particles read in total'
matter_data = np.asarray(matter_data,dtype=np.float32)
matter_file.close()
halo_file = open(haloDir+haloFilename,'r')
halo_data = np.reshape(halo_file.read().replace('\n',',').split(',')[0:-1],(-1,12)).astype(np.float32)
halo_file.close()
void_file = open(voidDir+voidFilename1,'r')
void_header1 = void_file.readline().split(",")
void_data1 = np.reshape(void_file.read().split(),(-1,len(void_header))).astype(np.float32)
void_file.close()
void_file = open(voidDir+voidFilename2,'r')
void_header2 = void_file.readline().split(",")+void_file.readline().split(",")
void_data2 = np.reshape(void_file.read().split(),(-1,11)).astype(np.float32)
void_file.close()
# Define arrays
xm = matter_data[:,0:3]
xh = halo_data[:,0:3]
mh = halo_data[:,6]
xv = void_data1[:,0:3]
rv = void_data1[:,4]
vv = void_data1[:,6]
mv = void_data2[:,8]*Mpart
# Interpolate to mesh # Interpolate to mesh
dm, wm, ws = xcorlib.cic(xm, Lbox, Lboxcut = Lboxcut, Nmesh = Nmesh) dm, wm, ws = xcorlib.cic(xm, Lbox, Lboxcut = Lboxcut, Nmesh = Nmesh, weights = None)
dh, wm, ws = xcorlib.cic(xh, Lbox, Lboxcut = Lboxcut, Nmesh = Nmesh) dv, wm, ws = xcorlib.cic(xv, Lbox, Lboxcut = Lboxcut, Nmesh = Nmesh, weights = None)
dv, wm, ws = xcorlib.cic(xv, Lbox, Lboxcut = Lboxcut, Nmesh = Nmesh)
# Load dark matter grid # Fourier transform
#output = open('dm_'+str(Nmesh)+'_ss'+str(ss)+'_z'+str(z)+'.dat', 'rb') dmk = np.fft.rfftn(dm)
#dm = pickle.load(output) dvk = np.fft.rfftn(dv)
#output.close()
# Save dark matter grid # 1D Power spectra & correlation functions
#output = open('dm_'+str(Nmesh)+'_ss'+str(ss)+'_z'+str(z)+'.dat', 'wb') ((Nm, km, Pmm, SPmm),(Nmx, rm, Xmm, SXmm)) = xcorlib.powcor(dmk, dmk, Lbox, Nbin, 'lin', True, True, 1)
#pickle.dump(dm,output) ((Nm, km, Pvm, SPvm),(Nmx, rm, Xvm, SXvm)) = xcorlib.powcor(dvk, dmk, Lbox, Nbin, 'lin', True, True, 1)
#output.close() ((Nm, km, Pvv, SPvv),(Nmx, rm, Xvv, SXvv)) = xcorlib.powcor(dvk, dvk, Lbox, Nbin, 'lin', True, True, 1)
# Power spectra & correlation functions
((Nm, km, Pmm, SPmm),(Nmx, rm, Xmm, SXmm)) = xcorlib.powcor(dm, dm, Lbox, Nmesh = Nmesh, Nbin = Nbin, scale = 'lin', cor = True)
((Nm, km, Pvm, SPvm),(Nmx, rm, Xvm, SXvm)) = xcorlib.powcor(dv, dm, Lbox, Nmesh = Nmesh, Nbin = Nbin, scale = 'lin', cor = True)
((Nm, km, Phm, SPhm),(Nmx, rm, Xhm, SXhm)) = xcorlib.powcor(dh, dm, Lbox, Nmesh = Nmesh, Nbin = Nbin, scale = 'lin', cor = True)
((Nm, km, Pvv, SPvv),(Nmx, rm, Xvv, SXvv)) = xcorlib.powcor(dv, dv, Lbox, Nmesh = Nmesh, Nbin = Nbin, scale = 'lin', cor = True)
((Nm, km, Pvh, SPvh),(Nmx, rm, Xvh, SXvh)) = xcorlib.powcor(dv, dh, Lbox, Nmesh = Nmesh, Nbin = Nbin, scale = 'lin', cor = True)
((Nm, km, Phh, SPhh),(Nmx, rm, Xhh, SXhh)) = xcorlib.powcor(dh, dh, Lbox, Nmesh = Nmesh, Nbin = Nbin, scale = 'lin', cor = True)
# 2D Power spectra & correlation functions
((Nm2d, kmper, kmpar, Pmm2d),(Nmx2d, rmper, rmpar, Xmm2d)) = xcorlib.powcor(dmk, dmk, Lbox, Nbin, 'lin', True, True, 2)
((Nm2d, kmper, kmpar, Pvm2d),(Nmx2d, rmper, rmpar, Xvm2d)) = xcorlib.powcor(dvk, dmk, Lbox, Nbin, 'lin', True, True, 2)
# Number densities # Number densities
nm = np.empty(len(km)) nm = np.empty(len(km))
nh = np.empty(len(km)) nh = np.empty(len(km))
nv = np.empty(len(km)) nv = np.empty(len(km))
nm[:] = Npart/Lbox**3 nm[:] = len(xm)/Lbox**3
nh[:] = len(xh)/Lbox**3 nh[:] = len(xh)/Lbox**3
nv[:] = len(xv)/Lbox**3 nv[:] = len(xv)/Lbox**3
# Number functions
Mbin = 40
Vbin = 40
Nh, Mh = np.histogram(mh, bins = mh.min()*(mh.max()/mh.min())**(np.arange(Mbin+1)/float(Mbin)))
Nvm, Mv = np.histogram(mv, bins = mv.min()*(mv.max()/mv.min())**(np.arange(Vbin+1)/float(Vbin)))
Nvv, Vv = np.histogram(vv, bins = vv.min()*(vv.max()/vv.min())**(np.arange(Vbin+1)/float(Vbin)))
# Bias # Bias
b_hh = np.sqrt(Phh/Pmm) b_hh = np.sqrt(Phh/Pmm)
b_vv = np.sqrt(Pvv/Pmm) b_vv = np.sqrt(Pvv/Pmm)
@ -189,109 +60,63 @@ def computeCrossCor(catalogDir,
b_vm = Pvm/Pmm b_vm = Pvm/Pmm
b_vh = Pvh/Phh b_vh = Pvh/Phh
knl = 0.04 # Wavenumber above which nonlinearities kick in [h/Mpc]
idls = np.where(km <= knl)[0]
bm_hm = np.average(b_hm[idls],weights=Nm[idls])
bm_vm = np.average(b_vm[idls],weights=Nm[idls])
bm_vh = np.average(b_vh[idls],weights=Nm[idls])
# Shot Noise # Shot Noise
sn_hh = Phh - Phm**2/Pmm sn_hh = Phh - Phm**2/Pmm
sn_vh = Pvh - Pvm*Phm/Pmm sn_vh = Pvh - Pvm*Phm/Pmm
sn_vv = Pvv - Pvm**2/Pmm sn_vv = Pvv - Pvm**2/Pmm
# Plots # Plots
ms = 4 mpl.rc('font', family='serif')
mew = 0.2 ms = 2.5
fs = 16
mew = 0.1
margin = 1.2 margin = 1.2
kmin = km.min()/margin kmin = km.min()/margin
kmax = km.max()*margin kmax = km.max()*margin
rmin = rm.min()/margin rmin = rm.min()/margin
rmax = rm.max()*margin rmax = rm.max()*margin
plt.imshow(np.sum(dm+1,2)/Nmesh,extent=[0,Lbox,0,Lbox],aspect='equal',cmap='YlGnBu_r',interpolation='gaussian')
# Density fields (projected)
plt.imshow(np.sum(dm[:,:,:]+1,2),extent=[0,Lbox,0,Lbox],aspect='equal',cmap='YlGnBu_r',interpolation='gaussian')
plt.xlabel(r'$x \;[h^{-1}\mathrm{Mpc}]$') plt.xlabel(r'$x \;[h^{-1}\mathrm{Mpc}]$')
plt.ylabel(r'$y \;[h^{-1}\mathrm{Mpc}]$') plt.ylabel(r'$y \;[h^{-1}\mathrm{Mpc}]$')
plt.title(r'Dark matter') plt.title(r'Dark matter')
plt.colorbar() plt.savefig(figDir+'/dm_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight")
plt.savefig(outputDir+'/dm_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight")
plt.clf() plt.clf()
plt.imshow(np.sum(dv+1,2)/Nmesh,extent=[0,Lbox,0,Lbox],aspect='equal',cmap='YlGnBu_r',interpolation='gaussian') plt.imshow(np.sum(dv[:,:,:]+1,2)/Nmesh,extent=[0,Lbox,0,Lbox],aspect='equal',cmap='YlGnBu_r',interpolation='gaussian')
plt.xlabel(r'$x \;[h^{-1}\mathrm{Mpc}]$') plt.xlabel(r'$x \;[h^{-1}\mathrm{Mpc}]$')
plt.ylabel(r'$y \;[h^{-1}\mathrm{Mpc}]$') plt.ylabel(r'$y \;[h^{-1}\mathrm{Mpc}]$')
plt.title(r'Voids') plt.title(r'Voids')
plt.colorbar() plt.savefig(figDir+'/dv_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight") #, dpi=300
plt.savefig(outputDir+'/dv_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight")
plt.clf()
plt.imshow(np.sum(dh+1,2)/Nmesh,extent=[0,Lbox,0,Lbox],aspect='equal',cmap='YlGnBu_r',interpolation='gaussian')
plt.xlabel(r'$x \;[h^{-1}\mathrm{Mpc}]$')
plt.ylabel(r'$y \;[h^{-1}\mathrm{Mpc}]$')
plt.title(r'Halos')
plt.colorbar()
plt.savefig(outputDir+'/dh_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight")
plt.clf() plt.clf()
pa, = plt.plot(Mh[:-1], Nh, 'ro-', ms=ms, mew=mew) # Power spectra & correlation functions
pb, = plt.plot(Vv[:-1]*1e9, Nvv, 'bo-', ms=ms, mew=mew) pa ,= plt.plot(km, Phh, 'r-s', ms=ms, mew=mew, mec='k')
plt.xlabel(r'$M \;[h^{-1}M_{\odot}]$ , $V \;[h^{-3}\mathrm{kpc}^3]$') plt.plot(km, Phh-sn_hh, 'r--', ms=ms, mew=mew)
plt.ylabel(r'$N(M,V)$') plt.plot(km, sn_hh, 'r:', ms=ms, mew=mew)
plt.title(r'Number of halos and voids')
plt.xscale('log')
plt.yscale('log')
plt.xlim(min(10**np.floor(np.log10(Vv.min()))*1e9,10**np.floor(np.log10(Mh.min()))), max(10**np.ceil(np.log10(Mh.max())),10**np.ceil(np.log10(Vv.max()))*1e8))
plt.ylim(10**np.floor(np.log10(Nh.min())), 10**np.ceil(np.log10(Nh.max())))
plt.annotate(r'$\frac{4\pi}{3}\langle r_\mathrm{v}\rangle^3$', xy=(4./3.*np.pi*rv.mean()**3*1e9,1.1), xytext=(-50,235),textcoords='offset points',arrowprops=dict(fc='k',arrowstyle="->",connectionstyle="angle,angleA=0,angleB=90,rad=10"))
plt.legend([pa,pb],['halos','voids'],'best' )
plt.savefig(outputDir+'/number_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight")
plt.clf()
plt.subplot(211)
plt.subplots_adjust(wspace=0,hspace=0)
plt.plot(np.sort(vv), mv[np.argsort(vv)], 'r-', ms=ms, mew=mew, lw=0.01)
plt.plot(np.sort(vv), np.sort(mv), 'k-', ms=ms, mew=mew)
plt.plot(np.sort(vv), np.sort(vv)*1e9, 'k--', ms=ms, mew=mew)
plt.title(r'Mass-volume relation of voids')
plt.ylabel(r'$M \;[h^{-1}M_\odot]$')
plt.xscale('log')
plt.yscale('log')
plt.subplot(211).xaxis.set_major_formatter(NullFormatter())
plt.subplot(212)
plt.plot(np.sort(vv), mv[np.argsort(vv)]/np.sort(vv)/rho_m-1., 'b-', ms=ms, mew=mew, lw=0.01)
plt.plot(np.sort(vv), np.sort(mv)/np.sort(vv)/rho_m-1., 'k-', ms=ms, mew=mew)
plt.xlabel(r'$V \;[h^{-3}\mathrm{Mpc}^3]$')
plt.ylabel(r'$\delta$')
plt.xscale('log')
plt.yscale('linear')
plt.ylim(-1.01,-0.861)
plt.savefig(outputDir+'/massvol_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight")
plt.clf()
pa ,= plt.plot(km, Phh, 'r-', ms=ms, mew=mew)
plt.plot(km, Phh-sn_hh, 'r:', ms=ms, mew=mew)
plt.fill_between(km, Phh+SPhh, abs(Phh-SPhh), color='r', alpha=0.2) plt.fill_between(km, Phh+SPhh, abs(Phh-SPhh), color='r', alpha=0.2)
pb ,= plt.plot(km, Phm, 'y-', ms=ms, mew=mew) pb ,= plt.plot(km, Phm, 'y-^', ms=ms, mew=mew, mec='k')
plt.fill_between(km, Phm+SPhm, abs(Phm-SPhm), color='y', alpha=0.2) plt.fill_between(km, Phm+SPhm, abs(Phm-SPhm), color='y', alpha=0.2)
pc ,= plt.plot(km, Pmm, 'k-', ms=ms, mew=mew) pc ,= plt.plot(km, Pmm, 'k-o', ms=0.8*ms, mew=mew, mec='k')
plt.plot(km, Pmm-1./nm, 'k:', ms=ms, mew=mew) plt.plot(km, Pmm-1./nm, 'k--', ms=ms, mew=mew)
plt.fill_between(km, Pmm+SPmm, abs(Pmm-SPmm), color='k', alpha=0.2) plt.fill_between(km, Pmm+SPmm, abs(Pmm-SPmm), color='k', alpha=0.2)
pd ,= plt.plot(km, Pvh, 'g-', ms=ms, mew=mew) pd ,= plt.plot(km, Pvh, 'g-*', ms=1.5*ms, mew=mew, mec='k')
plt.plot(km, -Pvh, 'g--', ms=ms, mew=mew) plt.plot(km, -Pvh, 'g*', ms=1.5*ms, mew=mew, mec='k')
plt.plot(km, abs(Pvh-sn_vh), 'g:', ms=ms, mew=mew) plt.plot(km, abs(Pvh-sn_vh), 'g--', ms=ms, mew=mew)
plt.plot(km, sn_vh, 'g:', ms=ms, mew=mew)
plt.plot(km, -sn_vh, 'g-.', ms=ms, mew=mew)
plt.fill_between(km, abs(Pvh+SPvh), abs(Pvh-SPvh), color='g', alpha=0.2) plt.fill_between(km, abs(Pvh+SPvh), abs(Pvh-SPvh), color='g', alpha=0.2)
pe ,= plt.plot(km, Pvm, 'm-', ms=ms, mew=mew) pe ,= plt.plot(km, Pvm, 'm-D', ms=ms, mew=mew, mec='k')
plt.plot(km, -Pvm, 'm--', ms=ms, mew=mew) plt.plot(km, -Pvm, 'mD', ms=ms, mew=mew, mec='k')
plt.fill_between(km, abs(Pvm+SPvm), abs(Pvm-SPvm), color='m', alpha=0.2) plt.fill_between(km, abs(Pvm+SPvm), abs(Pvm-SPvm), color='m', alpha=0.2)
pf ,= plt.plot(km, Pvv, 'b-', ms=ms, mew=mew) pf ,= plt.plot(km, Pvv, 'b-p', ms=1.3*ms, mew=mew, mec='k')
plt.plot(km, Pvv-sn_vv, 'b:', ms=ms, mew=mew) plt.plot(km, Pvv-sn_vv, 'b--', ms=ms, mew=mew)
plt.plot(km, sn_vv, 'b:', ms=ms, mew=mew)
plt.fill_between(km, Pvv+SPvv, abs(Pvv-SPvv), color='b', alpha=0.2) plt.fill_between(km, Pvv+SPvv, abs(Pvv-SPvv), color='b', alpha=0.2)
plt.annotate(r'$\frac{\pi}{\langle r_\mathrm{v}\rangle}$', xy=(np.pi/(rv.mean()),1.01*10**np.floor(np.log10(abs(Pvh).min()))/margin), xytext=(10,280),textcoords='offset points',arrowprops=dict(fc='k',arrowstyle="->",connectionstyle="angle,angleA=0,angleB=90,rad=10"))
plt.xlabel(r'$k \;[h\mathrm{Mpc}^{-1}]$') plt.xlabel(r'$k \;[h\mathrm{Mpc}^{-1}]$')
plt.ylabel(r'$P(k) \;[h^{-3}\mathrm{Mpc}^3]$') plt.ylabel(r'$P(k) \;[h^{-3}\mathrm{Mpc}^3]$')
plt.title(r'Power spectra') plt.title(r'Power spectra')
@ -299,11 +124,10 @@ def computeCrossCor(catalogDir,
plt.yscale('log') plt.yscale('log')
plt.xlim(kmin,kmax) plt.xlim(kmin,kmax)
plt.ylim(10**np.floor(np.log10(abs(Pvh).min()))/margin, max(10**np.ceil(np.log10(Phh.max())),10**np.ceil(np.log10(Pvv.max())))*margin) plt.ylim(10**np.floor(np.log10(abs(Pvh).min()))/margin, max(10**np.ceil(np.log10(Phh.max())),10**np.ceil(np.log10(Pvv.max())))*margin)
plt.legend([pa, pb, pc, pd, pe, pf],['hh', 'hm', 'mm', 'vh', 'vm', 'vv'],'lower left' ) plt.legend([pa, pb, pc, pd, pe, pf],['gg', 'gm', 'mm', 'vg', 'vm', 'vv'],'lower left',prop={'size':12})
plt.savefig(outputDir+'/power_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight") plt.savefig(figDir+'/power_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight")
plt.clf() plt.clf()
pa ,= plt.plot(rm, Xhh, 'r-', ms=ms, mew=mew) pa ,= plt.plot(rm, Xhh, 'r-', ms=ms, mew=mew)
plt.fill_between(rm, abs(Xhh+SXhh), abs(Xhh-SXhh), color='r', alpha=0.2) plt.fill_between(rm, abs(Xhh+SXhh), abs(Xhh-SXhh), color='r', alpha=0.2)
pb ,= plt.plot(rm, Xhm, 'y-', ms=ms, mew=mew) pb ,= plt.plot(rm, Xhm, 'y-', ms=ms, mew=mew)
@ -317,8 +141,8 @@ def computeCrossCor(catalogDir,
plt.plot(rm, -Xvm, 'm--', ms=ms, mew=mew) plt.plot(rm, -Xvm, 'm--', ms=ms, mew=mew)
plt.fill_between(rm, abs(Xvm+SXvm), abs(Xvm-SXvm), color='m', alpha=0.2) plt.fill_between(rm, abs(Xvm+SXvm), abs(Xvm-SXvm), color='m', alpha=0.2)
pf ,= plt.plot(rm, Xvv, 'b-', ms=ms, mew=mew) pf ,= plt.plot(rm, Xvv, 'b-', ms=ms, mew=mew)
plt.plot(rm, -Xvv, 'b--', ms=ms, mew=mew)
plt.fill_between(rm, abs(Xvv+SXvv), abs(Xvv-SXvv), color='b', alpha=0.2) plt.fill_between(rm, abs(Xvv+SXvv), abs(Xvv-SXvv), color='b', alpha=0.2)
plt.annotate(r'$\langle r_\mathrm{v}\rangle$', xy=(rv.mean(),1.01*10**np.floor(np.log10(abs(Xvh).min()))/margin), xytext=(10,300),textcoords='offset points',arrowprops=dict(fc='k',arrowstyle="->",connectionstyle="angle,angleA=0,angleB=90,rad=10"))
plt.xlabel(r'$r \;[h^{-1}\mathrm{Mpc}]$') plt.xlabel(r'$r \;[h^{-1}\mathrm{Mpc}]$')
plt.ylabel(r'$\xi(r)$') plt.ylabel(r'$\xi(r)$')
plt.title(r'Correlation functions') plt.title(r'Correlation functions')
@ -326,42 +150,51 @@ def computeCrossCor(catalogDir,
plt.yscale('log') plt.yscale('log')
plt.xlim(rmin,rmax) plt.xlim(rmin,rmax)
plt.ylim(10**np.floor(np.log10(abs(Xvh).min()))/margin, max(10**np.ceil(np.log10(Xhh.max())),10**np.ceil(np.log10(Xvv.max())))*margin) plt.ylim(10**np.floor(np.log10(abs(Xvh).min()))/margin, max(10**np.ceil(np.log10(Xhh.max())),10**np.ceil(np.log10(Xvv.max())))*margin)
plt.legend([pa, pb, pc, pd, pe, pf],['hh', 'hm', 'mm', 'vh', 'vm', 'vv'],'lower left' ) plt.legend([pa, pb, pc, pd, pe, pf],['gg', 'gm', 'mm', 'vg', 'vm', 'vv'],'best',prop={'size':12})
plt.savefig(outputDir+'/correlation_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight") plt.savefig(figDir+'/correlation_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight")
plt.clf() plt.clf()
pa, = plt.plot(km, b_hh, 'r-', ms=ms, mew=mew) # 2D power spectra & correlation functions
pb, = plt.plot(km, b_hm, 'r--', ms=ms, mew=mew) kpermin = kmper.min()
pc, = plt.plot(km, b_vv, 'b-', ms=ms, mew=mew) kpermax = 0.3001
pd, = plt.plot(km, b_vm, 'b--', ms=ms, mew=mew) kparmin = kmpar.min()
pe, = plt.plot(km, b_vh/bm_vh, 'g-', ms=ms, mew=mew) kparmax = 0.3001
plt.plot(km, np.sin(km*rv.mean())/(km*rv.mean()), 'k:', ms=ms, mew=mew) rpermin = rmper.min()
plt.xlabel(r'$k \;[h\mathrm{Mpc}^{-1}]$') rpermax = 40
plt.ylabel(r'$b(k)$') rparmin = rmpar.min()
plt.title(r'Bias') rparmax = 40
plt.xscale('log')
for (P2d,idx,vmin,vmax) in ([Pmm2d,'mm',None,None],[Pvm2d,'vm',None,None],[Phm2d,'gm',None,None],[Pvv2d,'vv',None,2.9],[Pvh2d,'vg',None,None],[Phh2d,'gg',None,None]):
cut = np.where(kmper[:,0] <= kpermax)[0].max()+2
plt.pcolormesh(kmper[0:cut,0:cut], kmpar[0:cut,0:cut], P2d[0:cut,0:cut]/1e4, cmap=cm.Spectral_r, shading='gouraud', vmin=vmin, vmax=vmax)
plt.colorbar(format='%.1f')
plt.contour(kmper[0:cut,0:cut], kmpar[0:cut,0:cut], P2d[0:cut,0:cut]/1e4, levels=np.array(P2d.min()+(P2d.max()-P2d.min())*(np.arange(Nbin/6+1)/float(Nbin/6)))/1e4, vmin=vmin, vmax=vmax, colors='k', linewidths=0.2)
plt.xlabel(r'$k_\perp \;[h\mathrm{Mpc}^{-1}]$', fontsize=fs)
plt.ylabel(r'$k_\parallel \;[h\mathrm{Mpc}^{-1}]$', fontsize=fs)
plt.axes().set_aspect('equal')
plt.xscale('linear')
plt.yscale('linear') plt.yscale('linear')
plt.xlim(kmin,kmax) plt.xlim(kpermin,kpermax)
plt.ylim(np.floor(b_vm.min()),np.ceil(max(b_hh.max(),b_vv.max()))) plt.ylim(kparmin,kparmax)
plt.legend([pa,pb,pc,pd,pe],['hh', 'hm', 'vv', 'vm', r'$\bar{u}_\mathrm{v}(k)$'],'best' ) plt.title(r'$P_{\mathrm{'+idx+r'}}(k_\perp, k_\parallel) \;[10^4h^{-3}\mathrm{Mpc}^3]$', fontsize=fs)
plt.savefig(outputDir+'/bias_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight") plt.savefig(figDir+'/P'+idx+'2d_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight")
plt.clf() plt.clf()
for (X2d,idx,vmin,vmax) in ([Xmm2d,'mm',None,None],[Xvm2d,'vm',None,None],[Xhm2d,'gm',None,None],[Xvv2d,'vv',None,0.2],[Xvh2d,'vg',None,None],[Xhh2d,'gg',None,None]):
pa, = plt.plot(km, sn_hh, 'r-', ms=ms, mew=mew) cut = np.where(rmper[:,0] <= rpermax)[0].max()+3
pb, = plt.plot(km, sn_vh, 'g-', ms=ms, mew=mew) plt.pcolormesh(rmper[0:cut,0:cut], rmpar[0:cut,0:cut], X2d[0:cut,0:cut], cmap=cm.Spectral_r, shading='gouraud', vmin=vmin, vmax=vmax)
pc, = plt.plot(km, sn_vv, 'b-', ms=ms, mew=mew) plt.colorbar(format='%+.2f')
plt.plot(km, abs(sn_vh), 'g:') plt.contour(rmper[0:cut,0:cut], rmpar[0:cut,0:cut], X2d[0:cut,0:cut], levels=np.array(X2d.min()+(X2d.max()-X2d.min())*(np.arange(Nbin/6+1)/float(Nbin/6))), vmin=vmin, vmax=vmax, colors='k', linewidths=0.2)
pd, = plt.plot(km, 1/nh, 'r--') plt.xlabel(r'$r_\perp \;[h^{-1}\mathrm{Mpc}]$', fontsize=fs)
pe, = plt.plot(km, 1/nv, 'b-.') plt.ylabel(r'$r_\parallel \;[h^{-1}\mathrm{Mpc}]$', fontsize=fs)
plt.xlabel(r'$k \;[h\mathrm{Mpc}^{-1}]$') plt.axes().set_aspect('equal')
plt.ylabel(r'$\sigma^2(k)$') plt.xscale('linear')
plt.title(r'Shotnoise') plt.yscale('linear')
plt.xscale('log') plt.xlim(rpermin,rpermax)
plt.yscale('log') plt.ylim(rparmin,rparmax)
plt.xlim(kmin,kmax) plt.title(r'$\xi_{\mathrm{'+idx+r'}}(r_\perp, r_\parallel)$', fontsize=fs)
plt.ylim(10**np.floor(np.log10(abs(sn_vh).min())), 10**np.ceil(np.log10(sn_vv.max()))) plt.savefig(figDir+'/X'+idx+'2d_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight")
plt.legend([pa,pb,pc,pd,pe],['hh', 'vh', 'vv', r'$\bar{n}_\mathrm{h}^{-1}$', r'$\bar{n}_\mathrm{v}^{-1}$'],'best' )
plt.savefig(outputDir+'/shotnoise_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight")
plt.clf() plt.clf()
return

View file

@ -0,0 +1,88 @@
import numpy as np
# CIC interpolation
def cic(x, Lbox, Lboxcut = 0, Nmesh = 128, weights = None):
if weights == None: weights = 1
wm = np.mean(weights)
ws = np.mean(weights**2)
d = np.mod(x/(Lbox+2*Lboxcut)*Nmesh,1)
box = ([Lboxcut,Lbox+Lboxcut],[Lboxcut,Lbox+Lboxcut],[Lboxcut,Lbox+Lboxcut])
rho = np.histogramdd(x, range = box, bins = Nmesh, weights = weights*(1-d[:,0])*(1-d[:,1])*(1-d[:,2]))[0] \
+ np.roll(np.histogramdd(x, range = box, bins = Nmesh, weights = weights*d[:,0]*(1-d[:,1])*(1-d[:,2]))[0],1,0) \
+ np.roll(np.histogramdd(x, range = box, bins = Nmesh, weights = weights*(1-d[:,0])*d[:,1]*(1-d[:,2]))[0],1,1) \
+ np.roll(np.histogramdd(x, range = box, bins = Nmesh, weights = weights*(1-d[:,0])*(1-d[:,1])*d[:,2])[0],1,2) \
+ np.roll(np.roll(np.histogramdd(x, range = box, bins = Nmesh, weights = weights*d[:,0]*d[:,1]*(1-d[:,2]))[0],1,0),1,1) \
+ np.roll(np.roll(np.histogramdd(x, range = box, bins = Nmesh, weights = weights*d[:,0]*(1-d[:,1])*d[:,2])[0],1,0),1,2) \
+ np.roll(np.roll(np.histogramdd(x, range = box, bins = Nmesh, weights = weights*(1-d[:,0])*d[:,1]*d[:,2])[0],1,1),1,2) \
+ np.roll(np.roll(np.roll(np.histogramdd(x, range = box, bins = Nmesh, weights = weights*d[:,0]*d[:,1]*d[:,2])[0],1,0),1,1),1,2)
rho /= wm
rho = rho/rho.mean() - 1.
return (rho, wm, ws)
# Power spectra & correlation functions
def powcor(d1, d2, Lbox, Nbin = 10, scale = 'lin', cor = False, cic = True, dim = 1):
Nmesh = len(d1)
# CIC correction
if cic:
wid = np.indices(np.shape(d1))
wid[np.where(wid >= Nmesh/2)] -= Nmesh
wid = wid*np.pi/Nmesh + 1e-100
wcic = np.prod(np.sin(wid)/wid,0)**2
# Shell average power spectrum
dk = 2*np.pi/Lbox
Pk = np.conj(d1)*d2*(Lbox/Nmesh**2)**3
if cic: Pk /= wcic**2
(Nm, km, Pkm, SPkm) = shellavg(np.real(Pk), dk, Nmesh, Nbin = Nbin, xmin = 0., xmax = Nmesh*dk/2, scale = scale, dim = dim)
# Inverse Fourier transform and shell average correlation function
if cor:
if cic: Pk *= wcic**2 # Undo cic-correction in correlation function
dx = Lbox/Nmesh
Xr = np.fft.irfftn(Pk)*(Nmesh/Lbox)**3
(Nmx, rm, Xrm, SXrm) = shellavg(np.real(Xr), dx, Nmesh, Nbin = Nbin/2, xmin = dx, xmax = 140., scale = scale, dim = dim)
return ((Nm, km, Pkm, SPkm),(Nmx, rm, Xrm, SXrm))
else: return (Nm, km, Pkm, SPkm)
# Shell averaging
def shellavg(f, dx, Nmesh, Nbin = 10, xmin = 0., xmax = 1., scale = 'lin', dim = 1):
x = np.indices(np.shape(f))
x[np.where(x >= Nmesh/2)] -= Nmesh
f = f.flatten()
if scale == 'lin': bins = xmin+(xmax-xmin)* np.linspace(0,1,Nbin+1)
if scale == 'log': bins = xmin*(xmax/xmin)**np.linspace(0,1,Nbin+1)
if dim == 1: # 1D
x = dx*np.sqrt(np.sum(x**2,0)).flatten()
Nm = np.histogram(x, bins = bins)[0]
xm = np.histogram(x, bins = bins, weights = x)[0]/Nm
fm = np.histogram(x, bins = bins, weights = f)[0]/Nm
fs = np.sqrt((np.histogram(x, bins = bins, weights = f**2)[0]/Nm - fm**2)/(Nm-1))
return (Nm, xm, fm, fs)
elif dim == 2: # 2D
xper = dx*np.sqrt(x[0,:,:,:]**2 + x[1,:,:,:]**2 + 1e-100).flatten()
xpar = dx*np.abs(x[2,:,:,:]).flatten()
x = dx*np.sqrt(np.sum(x**2,0)).flatten()
Nm = np.histogram2d(xper, xpar, bins = [bins,bins])[0]
xmper = np.histogram2d(xper, xpar, bins = [bins,bins], weights = xper)[0]/Nm
xmpar = np.histogram2d(xper, xpar, bins = [bins,bins], weights = xpar)[0]/Nm
fm = np.histogram2d(xper, xpar, bins = [bins,bins], weights = f)[0]/Nm
return (Nm, xmper, xmpar, fm)