diff --git a/python_tools/void_python_tools/voidUtil/xcorUtil.py b/python_tools/void_python_tools/voidUtil/xcorUtil.py
index 43f677f..9dac3c5 100644
--- a/python_tools/void_python_tools/voidUtil/xcorUtil.py
+++ b/python_tools/void_python_tools/voidUtil/xcorUtil.py
@@ -1,367 +1,200 @@
-#+
-#   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 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
-__all__=['computeCrossCor',]
-
-# ------------------------------------------------------------------------------
-
-def computeCrossCor(catalogDir, 
-                    outputDir="./", logDir="./",
-                    matchMethod="useID", dataPortion="central",
-                    strictMatch=True,
-                    pathToCTools="../../../c_tools"):
-
-# Computes void-void and void-matter(galaxy) correlations
-#  baseCatalogDir: directory of catalog 
-#  compareCatagDir: directory of catalog 2
-#  outputDir: directory to place outputs
-#  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):
-    os.makedirs(outputDir)
-
-  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
-  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()
+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
+import xcorlib
+
+def computeXcor(catalog,
+                figDir="./",
+                Nmesh = 256,
+                Nbin = 100
+                ):
+
+# Computes and plots void-void and void-matter(galaxy) correlations
+#   catalog: catalog to analyze
+#   figDir: where to place plots
+#   Nmesh: number of grid cells in power spectrum calculation
+#   Nbin: number of grid cells in plot
+
+  # Parameters
+  Lbox = catalog.boxLen # Boxlength
+  Lboxcut = 0.
+  Lbox -= 2*Lboxcut
+  
+  # Input particle arrays of shape (N,3)
+  xm = catalog.partPos # Halos / Galaxies / Dark matter
+  xv = catalog.voids[:].barycenter # Voids
+  
+  
+  # Interpolate to mesh
+  dm, wm, ws = xcorlib.cic(xm, Lbox, Lboxcut = Lboxcut, Nmesh = Nmesh, weights = None)
+  dv, wm, ws = xcorlib.cic(xv, Lbox, Lboxcut = Lboxcut, Nmesh = Nmesh, weights = None)
+  
+  # Fourier transform
+  dmk = np.fft.rfftn(dm)
+  dvk = np.fft.rfftn(dv)
+  
+  # 1D Power spectra & correlation functions
+  ((Nm, km, Pmm, SPmm),(Nmx, rm, Xmm, SXmm)) = xcorlib.powcor(dmk, dmk, Lbox, Nbin, 'lin', True, True, 1)
+  ((Nm, km, Pvm, SPvm),(Nmx, rm, Xvm, SXvm)) = xcorlib.powcor(dvk, dmk, Lbox, Nbin, 'lin', True, True, 1)
+  ((Nm, km, Pvv, SPvv),(Nmx, rm, Xvv, SXvv)) = xcorlib.powcor(dvk, dvk, Lbox, Nbin, 'lin', True, True, 1)
+  
+  # 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
+  nm = np.empty(len(km))
+  nh = np.empty(len(km))
+  nv = np.empty(len(km))
+  nm[:] = len(xm)/Lbox**3
+  nh[:] = len(xh)/Lbox**3
+  nv[:] = len(xv)/Lbox**3
+  
+  # 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
+  
+  # Shot Noise
+  sn_hh = Phh - Phm**2/Pmm
+  sn_vh = Pvh - Pvm*Phm/Pmm
+  sn_vv = Pvv - Pvm**2/Pmm
+  
+  
+  # Plots
+  mpl.rc('font', family='serif')
+  ms = 2.5
+  fs = 16
+  mew = 0.1
+  margin = 1.2
+  kmin = km.min()/margin
+  kmax = km.max()*margin
+  rmin = rm.min()/margin
+  rmax = rm.max()*margin
+  
+  
+  # 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.ylabel(r'$y \;[h^{-1}\mathrm{Mpc}]$')
+  plt.title(r'Dark matter')
+  plt.savefig(figDir+'/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.savefig(figDir+'/dv_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight") #, dpi=300
+  plt.clf()
+  
+  
+  # Power spectra & correlation functions
+  pa ,= plt.plot(km, Phh, 'r-s', ms=ms, mew=mew, mec='k')
+  plt.plot(km, Phh-sn_hh, 'r--', ms=ms, mew=mew)
+  plt.plot(km, 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, mec='k')
+  plt.fill_between(km, Phm+SPhm, abs(Phm-SPhm), color='y', alpha=0.2)
+  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.fill_between(km, Pmm+SPmm, abs(Pmm-SPmm), color='k', alpha=0.2)
+  pd ,= plt.plot(km, Pvh, 'g-*', ms=1.5*ms, mew=mew, mec='k')
+  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, 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)
+  pe ,= plt.plot(km, Pvm, 'm-D', ms=ms, mew=mew, mec='k')
+  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)
+  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, sn_vv, 'b:', ms=ms, mew=mew)
+  plt.fill_between(km, Pvv+SPvv, abs(Pvv-SPvv), color='b', alpha=0.2)
+  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],['gg', 'gm', 'mm', 'vg', 'vm', 'vv'],'lower left',prop={'size':12})
+  plt.savefig(figDir+'/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.plot(rm, -Xvv, 'b--', ms=ms, mew=mew)
+  plt.fill_between(rm, abs(Xvv+SXvv), abs(Xvv-SXvv), color='b', alpha=0.2)
+  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],['gg', 'gm', 'mm', 'vg', 'vm', 'vv'],'best',prop={'size':12})
+  plt.savefig(figDir+'/correlation_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight") 
+  plt.clf()
+  
+  
+  # 2D power spectra & correlation functions
+  kpermin = kmper.min()
+  kpermax = 0.3001
+  kparmin = kmpar.min()
+  kparmax = 0.3001
+  rpermin = rmper.min()
+  rpermax = 40
+  rparmin = rmpar.min()
+  rparmax = 40
+  
+  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.xlim(kpermin,kpermax)
+    plt.ylim(kparmin,kparmax)
+    plt.title(r'$P_{\mathrm{'+idx+r'}}(k_\perp, k_\parallel) \;[10^4h^{-3}\mathrm{Mpc}^3]$', fontsize=fs)
+    plt.savefig(figDir+'/P'+idx+'2d_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight")
+    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]):
+    cut = np.where(rmper[:,0] <= rpermax)[0].max()+3
+    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)
+    plt.colorbar(format='%+.2f')
+    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)  
+    plt.xlabel(r'$r_\perp \;[h^{-1}\mathrm{Mpc}]$', fontsize=fs)
+    plt.ylabel(r'$r_\parallel \;[h^{-1}\mathrm{Mpc}]$', fontsize=fs)
+    plt.axes().set_aspect('equal')
+    plt.xscale('linear')
+    plt.yscale('linear')
+    plt.xlim(rpermin,rpermax)
+    plt.ylim(rparmin,rparmax)
+    plt.title(r'$\xi_{\mathrm{'+idx+r'}}(r_\perp, r_\parallel)$', fontsize=fs)
+    plt.savefig(figDir+'/X'+idx+'2d_'+sample.fullName+'.pdf', format='pdf', bbox_inches="tight")
+    plt.clf()
+ 
+  return 
diff --git a/python_tools/void_python_tools/voidUtil/xcorlib.py b/python_tools/void_python_tools/voidUtil/xcorlib.py
new file mode 100644
index 0000000..3c370a8
--- /dev/null
+++ b/python_tools/void_python_tools/voidUtil/xcorlib.py
@@ -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)