> For the complete documentation index, see [llms.txt](https://docs.caterpillarproject.org/llms.txt). Markdown versions of documentation pages are available by appending `.md` to page URLs; this page is available as [Markdown](https://docs.caterpillarproject.org/usage/getting-started.md).

# Analyzing Data

## Analysis Modules

Two helpful libraries are our modules and analysis tools.

```bash
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:

```bash
setenv PYTHONPATH /path/to/modules:$PYTHONPATH
setenv PYTHONPATH /path/to/analysis:$PYTHONPATH
```

Lastly you will need to install `asciitable`, `h5py` and `pandas`.

{% hint style="success" %}
 These tools aren't critical, but they may make your life a lot easier.
{% endhint %}

## Halo Catalogs

Once you have the modules ready, uou can load a `ROCKSTAR` catalogue for a given snapshot simply as:

```python

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:

```python
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 Halo
  def 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 only
  def get_subs(self):
      # Get subhalo frame only
  def get_all_subs_recurse(self,haloID):
      # Retrieve all subhalos: sub and sub-sub, etc.
      # just need mask of all subhalos, then return data frame subset
  def get_all_subhalos_from_halo(self,haloID):
      # Retrieve all subhalos: sub and sub-sub, etc.
      # return pandas data frame of subhalos
  def get_all_sub_particles_from_halo(self,haloID):
      #returns int array of particle IDs belonging to all substructure
      #within host of haloID
  def get_all_particles_from_halo(self,haloID):
      #returns int array of all particles belonging to haloID
  def 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 7
  def 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_greg
  def H(self):
      #returns hubble parameter for rockstar run
  def get_most_gravbound_particles_from_halo(self,snapshot_dir, haloID):
      # Gets most bound particles just based on potential energy for specific halo ID
  def 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 halo
  def 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:

```python
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 id
parentid = 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 ids
zoomid = htils.load_zoomid(hpath)

# Return the pandas data frame with all the halos,
# return mvir for example
lastsnap = 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 masses
all_mvir = halos['mvir'] # units of Msol/h

# get the merger tree of the host and all its subs
cat = htils.load_mtc(hpath,haloids=[zoomid])

# read in every halo's merger tree
all_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 functions
mainbranch = 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:

```python
# load required modules
import haloutils as htils
import numpy as np

# select Cat-1 halo
hpaths = htils.get_paper_paths_lx(14)[0]

# select the last snapshot (z = 0)
snapshot = htils.get_lastsnap(hpath)

# load rockstar id of the host halo
zoomid = htils.load_zoomid(hpath)

#Load Halo Catalogue
halos = htils.load_rscat(hpath,snapshot)

#Select host halos
hosts = halos.get_hosts()

#Select subhalos
subs = halos.get_subs()

# Get positions of subs and hosts
print 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
```

## Merger Trees

Similarly the merger tree catalogues (once loaded) have a number of its own functions.

```python
def getMainBranch(self, row=0):
   """
   @param row: row of the halo you want the main branch for. Defaults to row 0
   Uses 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 massive
   These are all the subhalos that were destroyed in the interval
   """
```

These can be used in the following example:

```python
# load required modules

import haloutils as htils
import numpy as np

# select the caterpillar halo of interest
# based on halo id and resolution level

hid = 1387186
lx = 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 halo
zoomid = 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 catalogue
trees = 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 evolution
for mass,scale in zip(main_branch['mvir'],main_branch['scale']):
   print "%3.2f: %3.2e" % (scale,mass)

# output
1.00: 2.86e+14
0.99: 2.89e+14
0.98: 2.90e+14
0.97: 2.90e+14
0.95: 2.88e+14
0.94: 2.86e+14
0.93: 2.83e+14
0.92: 2.80e+14
0.91: 2.76e+14
0.90: 2.73e+14
0.89: 2.64e+14
0.88: 2.47e+14 ...

```

A function to find the descendant branch of any halo in merger tree catalogue. You should use it as follows:

```python
# get a tree of interest
mtc = haloutils.load_zoom_mtc(hpath)
host = mtc.Trees[0]

# make a dictionary that maps ids to rows
desc_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)
```

## Particle Data

If you want the Gadget header:

```python
# get the first halo in the catalogue

hpaths = htils.get_paper_paths_lx(14)[0]

# read its Gadget header
header = htils.get_halo_header(hpath)

# header contains the typical Gadget info
header.boxsize        header.massarr        header.omegaL
header.cooling        header.metals         header.redshift
header.double         header.nall           header.sfr
header.feedback       header.nall_highword  header.stellar_age
header.filenum        header.npart          header.time
header.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).

```python
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:

```python
import haloutils as htils

hid = 1387186
lx = 14
hpath = htils.hid_hpath_lx(hid,lx)

pos = htils.load_partblock(hpath,319,"POS ")
mass = htils.load_partblock(hpath,319,"MASS")
```

{% hint style="info" %}
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.
{% endhint %}

If you wanted just the ids for a selection of particle ids:

```python
pos = htils.load_partblock(hpath,zoomid,"POS ",partids=[listofids]) # units Mpc/h
```