Two helpful libraries are our modules and analysis tools.
git clone git@github.com:caterpillarproject/modules.git # Python 2.7+git clone git@github.com:caterpillarproject/analysis.git # Python 2.7+
Add these to your PYTHONPATH
environment variable, e.g. for .cshrc
add:
setenv PYTHONPATH /path/to/modules:$PYTHONPATHsetenv PYTHONPATH /path/to/analysis:$PYTHONPATH
Lastly you will need to install asciitable
, h5py
and pandas
.
These tools aren't critical, but they may make your life a lot easier.
Once you have the modules ready, uou can load a ROCKSTAR
catalogue for a given snapshot simply as:
​rscat = htils.load_rscat(hpath,319,verbose=True) # snapshot = 319 (z = 0)
Once you do this however, you will have access to the following methods:
def __init__(self, dir, snap_num, version=2, sort_by='mvir', base='halos_', digits=2, AllParticles=False):def get_particles_from_halo(self, haloID):# @param haloID: id number of halo. Not its row position in matrix# @return: a list of particle IDs in the Halodef get_subhalos_from_halo(self,haloID):#Retrieve subhalos only one level deep.#Does not get sub-sub halos, etc.def get_subhalos_from_halos(self,haloIDs):#Returns an array of pandas data frames of subhalos. one data frame#for each host halo. returns only first level of subhalos.def get_subhalos_from_halos_flat(self,haloIDs):#Returns a flattened pandas data frame of all subhalos within#the hosts given by haloIDs. Returns only first level of subhalos.def get_hosts(self):# Get host halo frame onlydef get_subs(self):# Get subhalo frame onlydef get_all_subs_recurse(self,haloID):# Retrieve all subhalos: sub and sub-sub, etc.# just need mask of all subhalos, then return data frame subsetdef get_all_subhalos_from_halo(self,haloID):# Retrieve all subhalos: sub and sub-sub, etc.# return pandas data frame of subhalosdef get_all_sub_particles_from_halo(self,haloID):#returns int array of particle IDs belonging to all substructure#within host of haloIDdef get_all_particles_from_halo(self,haloID):#returns int array of all particles belonging to haloIDdef get_all_num_particles_from_halo(self,haloID):# Get the actual number of particles 'total_npart' from halo as opposed to 'npart'.# mainly for versions less than 7def get_block_from_halo(self, snapshot_dir, haloID, blockname, allparticles=True):# quick load a block (hdf5 block) of particles belong to halo.# e.g. you want particle positions for haloid = 10 (use blockname="pos")# this works fastest on snapshots ordered by id and requires import readsnapHDF5_gregdef H(self):#returns hubble parameter for rockstar rundef get_most_gravbound_particles_from_halo(self,snapshot_dir, haloID):# Gets most bound particles just based on potential energy for specific halo IDdef get_most_bound_particles_from_halo(self, snapshot_dir, haloID):# Gets most bound particles for halo based on pot. energy and kin. energy# if potential block does not exist, it is calculate assuming a spherical halodef getversion(self):# returns the version of rockstar the run was done within# this will include versions made by Alex Ji, Greg Dooley & Brendan Griffen
One workflow might look as follows:
import haloutils as htils​# load the first 24 halos# (just change 14 to 11 for the lower resolution halos)​hpaths = htils.get_paper_paths_lx(14)hpath = hpaths[0] # select first Caterpillar halo​# strip down the path to just the halo idparentid = htils.get_parent_hid(hpath)​# get the central host id of the zoom-in halo# parentid is named from the parent simulation# the zooms have different idszoomid = htils.load_zoomid(hpath)​# Return the pandas data frame with all the halos,# return mvir for examplelastsnap = htils.get_lastsnap(hpath)halos = htils.load_rscat(hpath,lastsnap,verbose=True)mvir_host = halos.ix[zoomid]['mvir'] # units of Msol/h​# get one specific parameter of the host at# z = 0 (quick version of above)mvir_host = htils.get_quant_zoom(hpath,'mvir') # units of Msol/h​# return all halo virial massesall_mvir = halos['mvir'] # units of Msol/h​# get the merger tree of the host and all its subscat = htils.load_mtc(hpath,haloids=[zoomid])​# read in every halo's merger treeall_trees = htils.load_mtc(hpath,indexbyrsid=True)tree = all_trees[zoomid]​# if you feed more than one id to the above# it would be cat[0],cat[1] etc.​# you can now access the main branch > try tree. then tab complete to see other functionsmainbranch = tree.getMainBranch()​# we also have a short hand version:# mainbranch = htils.get_mainbranch(hpath)# which is very fast and skips the reading of the entire progenitor tree if you don't need dit# see further down this page for more options# output the main branch virial mass​print mainbranch['mvir']
Here is an example of some of these in action:
# load required modulesimport haloutils as htilsimport numpy as np​# select Cat-1 halohpaths = htils.get_paper_paths_lx(14)[0]​# select the last snapshot (z = 0)snapshot = htils.get_lastsnap(hpath)​# load rockstar id of the host halozoomid = htils.load_zoomid(hpath)​#Load Halo Cataloguehalos = htils.load_rscat(hpath,snapshot)​#Select host haloshosts = halos.get_hosts()​#Select subhalossubs = halos.get_subs()​# Get positions of subs and hostsprint hosts[['posX','posY','posZ']]print subs[['posX','posY','posZ']]​#Get particle ids from halo of interest (here it is the host)print halos.get_particles_from_halo(zoomid)​#Get virial radius of a specific halo id (in this case the host)print halos.ix[zoomid]['rvir'] # units of kpc/h
Similarly the merger tree catalogues (once loaded) have a number of its own functions.
def getMainBranch(self, row=0):"""@param row: row of the halo you want the main branch for. Defaults to row 0Uses getSubTree, then finds the smallest dfid that has no progenitors@return: all halos that are in the main branch of the halo specified by row (in a np structured array)"""def getMMP(self, row):"""@param row: row number (int) of halo considered@return: row number of the most massive parent, or None if no parent"""def getNonMMPprogenitors(self,row):"""return row index of all progenitors that are not the most massiveThese are all the subhalos that were destroyed in the interval"""
These can be used in the following example:
# load required modules​import haloutils as htilsimport numpy as np​# select the caterpillar halo of interest# based on halo id and resolution level​hid = 1387186lx = 14​hpath = htils.hid_hpath_lx(hid,lx)​# select the last snapshot (z = 0)snapshot = htils.get_lastsnap(hpath)​# load rockstar id of the host halozoomid = htils.load_zoomid(hpath)​# Load every tree (VERY slow)trees = htils.load_mtc(hpath)​# Just look at the first tree# tree = cat[0]# You can access all trees via: cat[0], cat[1] etc.​# If you want to load the tree for a particular ID from the rockstar cataloguetrees = htils.load_mtc(haloids=[zoomid]) # in this case the host (quite slow)tree = cat[0] # tree will contain all progenitors, including subhalos​# Just say you want to index the merger tree by the z = 0 root rockstar id# (i.e. the base of the tree). This is quite powerful because you might select# halos of interest in the rockstar catalogue then want to know, just for those# what their merger tree is (e.g. say you just want dwarf systems of a particular size)​trees = htils.load_mtc(haloids=[zoomid],indexbyrsid=True)tree = trees[zoomid]​# You can just loop through rockstar ids (if you gave it more than one id above)# and get out the accretion histories for a small sample of trees quite quickly​# to get the main branch​main_branch = tree.getMainBranch()​# print mass evolutionfor mass,scale in zip(main_branch['mvir'],main_branch['scale']):print "%3.2f: %3.2e" % (scale,mass)​# output1.00: 2.86e+140.99: 2.89e+140.98: 2.90e+140.97: 2.90e+140.95: 2.88e+140.94: 2.86e+140.93: 2.83e+140.92: 2.80e+140.91: 2.76e+140.90: 2.73e+140.89: 2.64e+140.88: 2.47e+14 ...​
A function to find the descendant branch of any halo in merger tree catalogue. You should use it as follows:
# get a tree of interestmtc = haloutils.load_zoom_mtc(hpath)host = mtc.Trees[0]​# make a dictionary that maps ids to rowsdesc_map = host.get_desc_map()​# get the descendent branch.desc_branch = host.getDescBranch(row, desc_map)​# similarly for main branches, you can use a dictionary that speeds up the# getMainBranch call substantially.mmp_map = host.get_mmp_map()main_branch = host.getMainBranch(row, mmp_map)​# you can get the branches without making the map as below, but they will be much slower.# only faster if you use ~ < 10 calls to the function.desc_branch = host.getDescBranch(row)main_branch = host.getMainBranch(row)
If you want the Gadget header:
# get the first halo in the catalogue​hpaths = htils.get_paper_paths_lx(14)[0]​# read its Gadget headerheader = htils.get_halo_header(hpath)​# header contains the typical Gadget infoheader.boxsize header.massarr header.omegaLheader.cooling header.metals header.redshiftheader.double header.nall header.sfrheader.feedback header.nall_highword header.stellar_ageheader.filenum header.npart header.timeheader.hubble header.omega0
Be sure to divide the relevant quantities (pos, rvir etc.) by header.hubble
. See the Gadget section in the sidebar for more information on the header and block types available.
If you wanted to get the postions of all the particles for a specific halo (or any block).
pos = htils.load_partblock(hpath,zoomid,"POS ") # units Mpc/h# "VEL ", "ID ", "MASS" etc. also work (notice the space)# check readsnapshots/readsnapHDF5.py for the other block names# you can call in the caterpillar modules​print pos*1000. # kpc/h​# output[ 32085.45117188 57312.79296875 44314.35546875][ 27002.18554688 10062.73242188 9899.70019531][ 26711.08789062 9560.22460938 10165.18847656][ 49757.3515625 21470.00195312 6461.90917969]...
If you want to read in the entire block, use the following:
import haloutils as htils​hid = 1387186lx = 14hpath = htils.hid_hpath_lx(hid,lx)​pos = htils.load_partblock(hpath,319,"POS ")mass = htils.load_partblock(hpath,319,"MASS")
Note that the mass block will have different values depending on how many layers of refinement there are for that zoom in simulation. If you use this code on a parent simulation it will be an array of length N all of the same value because there is only one particle type.
If you wanted just the ids for a selection of particle ids:
pos = htils.load_partblock(hpath,zoomid,"POS ",partids=[listofids]) # units Mpc/h