reading in Websky's halos.pksc

[borisbolliet]

Hi,

I think the file nemo/examples/SOSims/readWebSkyInputCatalog.py is not synchronized with latest Nemo's commits (makeACTName not found), and also I had to copy paste the web sky/cosmology.py file to get zofchi.
Attached is the modified file that runs for me:

"""

Reads in the halos.pksc file and outputs .fits and .reg format catalogs for
the 1 million most massive halos (M200m > 1e14 MSun). Cuts down to SO dec
range.

"""
import os
import sys
import healpy as hp
import astropy.table as atpy
from nemo import catalogs
# from cosmology import *
import numpy as np
import IPython


from scipy.interpolate import *
import numpy as np

omegab = 0.049
omegac = 0.261
omegam = omegab + omegac
h      = 0.68
ns     = 0.965
sigma8 = 0.81

c = 3e5

H0 = 100*h
nz = 100000
z1 = 0.0
z2 = 6.0
za = np.linspace(z1,z2,nz)
dz = za[1]-za[0]

H      = lambda z: H0*np.sqrt(omegam*(1+z)**3+1-omegam)
dchidz = lambda z: c/H(z)

chia = np.cumsum(dchidz(za))*dz

zofchi = interp1d(chia,za)


def makeACTName(RADeg, decDeg, prefix = 'ACT-CL'):
    """Makes ACT cluster name from RADeg, decDeg

    """

    actName=prefix+" J"+makeRA(RADeg)+makeDec(decDeg)

    return actName

def makeRA(myRADeg):
    """Makes RA part of ACT names.

    """
    hours=(myRADeg/360)*24
    if hours<10:
        sHours="0"+str(hours)[0]
    else:
        sHours=str(hours)[:2]

    mins=float(str(hours)[str(hours).index("."):])*60
    if mins<10:
        sMins="0"+str(mins)[:3]
    else:
        sMins=str(mins)[:4]

    return (sHours+sMins)#[:-2] # Trims off .x as not used in ACT names

#------------------------------------------------------------------------------------------------------------
def makeDec(myDecDeg):
    """Makes dec part of ACT names

    """

    # Positive
    if myDecDeg>0:
        if myDecDeg<10:
            sDeg="0"+str(myDecDeg)[0]
        else:
            sDeg=str(myDecDeg)[:2]

        mins=float(str(myDecDeg)[str(myDecDeg).index("."):])*60
        if mins<10:
            sMins="0"+str(mins)[:1]
        else:
            sMins=str(mins)[:2]

        return "+"+sDeg+sMins
    else:
        if myDecDeg>-10:
            sDeg="-0"+str(myDecDeg)[1]
        else:
            sDeg=str(myDecDeg)[:3]

        mins=float(str(myDecDeg)[str(myDecDeg).index("."):])*60
        if mins<10:
            sMins="0"+str(mins)[:1]
        else:
            sMins=str(mins)[:2]

        return str(sDeg+sMins)



rho = 2.775e11*omegam*h**2 # Msun/Mpc^3

# f=open('../WebSky/halos.pksc', 'rb')
f=open('/path/to/halos.pksc','rb')
N=np.fromfile(f,count=3,dtype=np.int32)[0]

# only take first five entries for testing (there are ~8e8 halos total...)
# comment the following line to read in all halos
N = 1000000

catalog=np.fromfile(f,count=N*10,dtype=np.float32)
catalog=np.reshape(catalog,(N,10))

x  = catalog[:,0];  y = catalog[:,1];  z = catalog[:,2] # Mpc (comoving)
vx = catalog[:,3]; vy = catalog[:,4]; vz = catalog[:,5] # km/sec
R  = catalog[:,6] # Mpc

# convert to mass, comoving distance, radial velocity, redshfit, RA and DEc
M        = 4*np.pi/3.*rho*R**3        # Msun (M200m)
chi      = np.sqrt(x**2+y**2+z**2)    # Mpc
vrad     = (x*vx + y*vy + z*vz) / chi # km/sec
redshift = zofchi(chi)

theta, phi  = hp.vec2ang(np.column_stack((x,y,z))) # in radians

decDeg=-1*(np.degrees(theta)-90) # Because HEALPix
RADeg=np.degrees(phi)

names=[]
for ra, dec in zip(RADeg, decDeg):
    names.append(makeACTName(ra, dec, prefix= 'MOCK-CL'))

outFileName="halos.fits"
tab=atpy.Table()
tab.add_column(atpy.Column(names, 'name'))
tab.add_column(atpy.Column(RADeg, 'RADeg'))
tab.add_column(atpy.Column(decDeg, 'decDeg'))
tab.add_column(atpy.Column(M, 'M200m'))
tab.add_column(atpy.Column(redshift, "z"))
tab=tab[np.where(tab['decDeg'] < 30)]
tab=tab[np.where(tab['decDeg'] > -70)]
tab.write(outFileName, overwrite = True)

catalogs.catalog2DS9(tab, "halos.reg")

#IPython.embed()
#sys.exit()

### e.g. project to a map, matching the websky orientations
#nside = 1024
#map   = np.zeros((hp.nside2npix(nside)))

#pix = hp.vec2pix(nside, x, y, z)
#pix = hp.ang2pix(nside, theta, phi) does the same

#weight = 1. #1 for number density, array of size(x) for arbitrary
#np.add.at(map, pix, weight)