From f699ed372f6791fabf250454a45ec974774f59f8 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Fri, 22 Mar 2013 11:46:20 -0500 Subject: [PATCH] added nico's cross-correlation script to a new directory, analysis/ --- analysis/datasetsToAnalyze.py | 51 +++ analysis/xcor.py | 352 ++++++++++++++++++ pipeline/generateCatalog.py | 2 +- .../void_python_tools/plotting/plotTools.py | 2 +- 4 files changed, 405 insertions(+), 2 deletions(-) create mode 100644 analysis/datasetsToAnalyze.py create mode 100755 analysis/xcor.py diff --git a/analysis/datasetsToAnalyze.py b/analysis/datasetsToAnalyze.py new file mode 100644 index 0000000..412de5c --- /dev/null +++ b/analysis/datasetsToAnalyze.py @@ -0,0 +1,51 @@ +#+ +# VIDE -- Void IDEntification pipeline -- ./crossCompare/analysis/datasetsToAnalyze.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. +#+ +#!/usr/bin/env python + + +outputDir = "/home/psutter2/workspace/Voids/analysis/xcor/" + +# Sim parameters +Ni = 1000237 # Number of dark matter particles per dimension in simulation +ss = 0.1 # Subsampling fraction of dark matter particles to read +Mpart = 8.721e9 # Particle mass [M_sol] +Npart = int(ss*Ni) # Particle number of subsample +Lboxcut = 0. # Size of optional margin to be cut from the box [h^(-1)Mpc] +Nmesh = 256 # Interpolation meshlength +Nbin = 70 # Number of bins for power spectrum and correlation function +r_H = 3000. # Hubble scale [h^(-1)Mpc] +ns = 0.95 # Spectral index +sigma_8 = 0.82 # Sigma_8 +h = 0.7 # Dimensionless Hubble parameter + +# Input files +matterDir = '/home/psutter2/workspace/Voids/catalogs/mergertree512/' +haloDir = '/home/psutter2/workspace/Voids/catalogs/mergertree512/' +matterFilename = 'mf_4s_1G_512_1.000' +haloFilename = 'mf_4s_1G_512_bgc2_1.000.sdf' +voidBaseDir = "/home/psutter2/workspace/Voids/" + + +sampleDirList = [ + "mergertree512/mt_ss0.01/sample_mt_ss0.01_z0.00_d00/", + ] + +dataPortion = "central" + diff --git a/analysis/xcor.py b/analysis/xcor.py new file mode 100755 index 0000000..2343443 --- /dev/null +++ b/analysis/xcor.py @@ -0,0 +1,352 @@ +#!/usr/bin/env python +#+ +# VIDE -- Void IDEntification pipeline -- ./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 * +from void_python_tools.plotting import * +import void_python_tools.xcor as xcorlib +import imp +import pickle +import argparse +import os +import string +import numpy as np +import matplotlib as mpl +mpl.use('Agg') +import matplotlib.pyplot as plt +import matplotlib.cm as cm +from matplotlib import rc +from matplotlib.ticker import NullFormatter +import random +import sys + +# ------------------------------------------------------------------------------ + +dataNameBase = "xcor" + +parser = argparse.ArgumentParser(description='Analyze.') +parser.add_argument('--parmFile', dest='parmFile', default='datasetsToAnalyze.py', + help='path to parameter file') +args = parser.parse_args() + +# ------------------------------------------------------------------------------ + +filename = args.parmFile +print " Loading parameters from", filename +if not os.access(filename, os.F_OK): + print " Cannot find parameter file %s!" % filename + exit(-1) +parms = imp.load_source("name", filename) +globals().update(vars(parms)) + +if not os.access(outputDir, os.F_OK): + os.makedirs(outputDir) + +for (iSample, sampleDir) in enumerate(sampleDirList): + + with open(voidBaseDir+sampleDir+"/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 + dm, wm, ws = xcorlib.cic(xm, Lbox, Lboxcut = Lboxcut, Nmesh = Nmesh) + dh, wm, ws = xcorlib.cic(xh, Lbox, Lboxcut = Lboxcut, Nmesh = Nmesh) + dv, wm, ws = xcorlib.cic(xv, Lbox, Lboxcut = Lboxcut, Nmesh = Nmesh) + + # Load dark matter grid + #output = open('dm_'+str(Nmesh)+'_ss'+str(ss)+'_z'+str(z)+'.dat', 'rb') + #dm = pickle.load(output) + #output.close() + + # Save dark matter grid + #output = open('dm_'+str(Nmesh)+'_ss'+str(ss)+'_z'+str(z)+'.dat', 'wb') + #pickle.dump(dm,output) + #output.close() + + + # 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) + + + # Number densities + nm = np.empty(len(km)) + nh = np.empty(len(km)) + nv = np.empty(len(km)) + nm[:] = Npart/Lbox**3 + nh[:] = len(xh)/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 + b_hh = np.sqrt(Phh/Pmm) + b_vv = np.sqrt(Pvv/Pmm) + b_hm = Phm/Pmm + b_vm = Pvm/Pmm + 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 + sn_hh = Phh - Phm**2/Pmm + sn_vh = Pvh - Pvm*Phm/Pmm + sn_vv = Pvv - Pvm**2/Pmm + + + + # Plots + ms = 4 + mew = 0.2 + margin = 1.2 + kmin = km.min()/margin + kmax = km.max()*margin + rmin = rm.min()/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') + plt.xlabel(r'$x \;[h^{-1}\mathrm{Mpc}]$') + plt.ylabel(r'$y \;[h^{-1}\mathrm{Mpc}]$') + plt.title(r'Dark matter') + plt.colorbar() + plt.savefig(outputDir+'/dm_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight") + plt.clf() + + 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.ylabel(r'$y \;[h^{-1}\mathrm{Mpc}]$') + plt.title(r'Voids') + plt.colorbar() + 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() + + + pa, = plt.plot(Mh[:-1], Nh, 'ro-', ms=ms, mew=mew) + pb, = plt.plot(Vv[:-1]*1e9, Nvv, 'bo-', ms=ms, mew=mew) + plt.xlabel(r'$M \;[h^{-1}M_{\odot}]$ , $V \;[h^{-3}\mathrm{kpc}^3]$') + plt.ylabel(r'$N(M,V)$') + 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) + pb ,= plt.plot(km, Phm, 'y-', ms=ms, mew=mew) + plt.fill_between(km, Phm+SPhm, abs(Phm-SPhm), color='y', alpha=0.2) + pc ,= plt.plot(km, Pmm, '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) + pd ,= plt.plot(km, Pvh, 'g-', ms=ms, mew=mew) + plt.plot(km, -Pvh, 'g--', ms=ms, mew=mew) + plt.plot(km, abs(Pvh-sn_vh), 'g:', ms=ms, mew=mew) + 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) + plt.plot(km, -Pvm, 'm--', ms=ms, mew=mew) + 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) + plt.plot(km, Pvv-sn_vv, 'b:', ms=ms, mew=mew) + 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.ylabel(r'$P(k) \;[h^{-3}\mathrm{Mpc}^3]$') + plt.title(r'Power spectra') + plt.xscale('log') + plt.yscale('log') + 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.legend([pa, pb, pc, pd, pe, pf],['hh', 'hm', 'mm', 'vh', 'vm', 'vv'],'lower left' ) + plt.savefig(outputDir+'/power_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight") + plt.clf() + + + 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) + pb ,= plt.plot(rm, Xhm, 'y-', ms=ms, mew=mew) + plt.fill_between(rm, abs(Xhm+SXhm), abs(Xhm-SXhm), color='y', alpha=0.2) + pc ,= plt.plot(rm, Xmm, 'k-', ms=ms, mew=mew) + plt.fill_between(rm, abs(Xmm+SXmm), abs(Xmm-SXmm), color='k', alpha=0.2) + pd ,= plt.plot(rm, Xvh, 'g-', ms=ms, mew=mew) + plt.plot(rm, -Xvh, 'g--', ms=ms, mew=mew) + plt.fill_between(rm, abs(Xvh+SXvh), abs(Xvh-SXvh), color='g', alpha=0.2) + pe ,= 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) + pf ,= 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.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.ylabel(r'$\xi(r)$') + plt.title(r'Correlation functions') + plt.xscale('log') + plt.yscale('log') + 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.legend([pa, pb, pc, pd, pe, pf],['hh', 'hm', 'mm', 'vh', 'vm', 'vv'],'lower left' ) + plt.savefig(outputDir+'/correlation_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight") + plt.clf() + + + pa, = plt.plot(km, b_hh, 'r-', ms=ms, mew=mew) + pb, = plt.plot(km, b_hm, 'r--', ms=ms, mew=mew) + pc, = plt.plot(km, b_vv, 'b-', ms=ms, mew=mew) + pd, = plt.plot(km, b_vm, 'b--', ms=ms, mew=mew) + pe, = plt.plot(km, b_vh/bm_vh, 'g-', ms=ms, mew=mew) + plt.plot(km, np.sin(km*rv.mean())/(km*rv.mean()), 'k:', ms=ms, mew=mew) + plt.xlabel(r'$k \;[h\mathrm{Mpc}^{-1}]$') + plt.ylabel(r'$b(k)$') + plt.title(r'Bias') + plt.xscale('log') + plt.yscale('linear') + plt.xlim(kmin,kmax) + plt.ylim(np.floor(b_vm.min()),np.ceil(max(b_hh.max(),b_vv.max()))) + plt.legend([pa,pb,pc,pd,pe],['hh', 'hm', 'vv', 'vm', r'$\bar{u}_\mathrm{v}(k)$'],'best' ) + plt.savefig(outputDir+'/bias_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight") + plt.clf() + + + pa, = plt.plot(km, sn_hh, 'r-', ms=ms, mew=mew) + pb, = plt.plot(km, sn_vh, 'g-', ms=ms, mew=mew) + pc, = plt.plot(km, sn_vv, 'b-', ms=ms, mew=mew) + plt.plot(km, abs(sn_vh), 'g:') + pd, = plt.plot(km, 1/nh, 'r--') + pe, = plt.plot(km, 1/nv, 'b-.') + plt.xlabel(r'$k \;[h\mathrm{Mpc}^{-1}]$') + plt.ylabel(r'$\sigma^2(k)$') + plt.title(r'Shotnoise') + plt.xscale('log') + plt.yscale('log') + plt.xlim(kmin,kmax) + plt.ylim(10**np.floor(np.log10(abs(sn_vh).min())), 10**np.ceil(np.log10(sn_vv.max()))) + 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() diff --git a/pipeline/generateCatalog.py b/pipeline/generateCatalog.py index 5ef2771..937a628 100755 --- a/pipeline/generateCatalog.py +++ b/pipeline/generateCatalog.py @@ -29,7 +29,7 @@ import pickle # ------------------------------------------------------------------------------ if (len(sys.argv) == 1): - print "Usage: ./analyzeVoids.py parameter_file.py" + print "Usage: ./generateCatalog.py parameter_file.py" exit(-1) if (len(sys.argv) > 1): diff --git a/python_tools/void_python_tools/plotting/plotTools.py b/python_tools/void_python_tools/plotting/plotTools.py index 317efef..b7d3efa 100644 --- a/python_tools/void_python_tools/plotting/plotTools.py +++ b/python_tools/void_python_tools/plotting/plotTools.py @@ -305,7 +305,7 @@ def plotVoidDistribution(workDir=None, sampleList=None, figDir=None, for (iSample,sample) in enumerate(sampleList): plt.clf() plt.xlabel("x (Mpc/h)") - plt.xlabel("y (Mpc/h)") + plt.ylabel("y (Mpc/h)") sampleName = sample.fullName