Examples

The examples shown here, except for the last one, use built-in functions found in readtreeHDF5.py, such as: The names of these functions attempt to be self-explanatory, but brief descriptions can be found inside the code.

The Python script for reading trees in "database mode" can be found here:

/n/home10/vrodrigu/Python/PythonModules/OwnModules/readtreeHDF5.py

Example 1

Print the snapshot number, stellar mass, and star formation rate along the main branch of a subhalo. Note that the keysel argument is used to specify which fields from the Subfind catalogs should be loaded. 
import readtreeHDF5
treedir = '/n/ghernquist/vrodrigu/MergerTrees/output/Subhalos/Illustris/L75n1820FP'
tree = readtreeHDF5.TreeDB(treedir)
snapnum = 135; subfind_id = 1
branch = tree.get_main_branch(snapnum, subfind_id, keysel=['SnapNum', 'SubhaloMassType', 'SubhaloSFR'])
for i in range(branch.nrows):
    print branch.SnapNum[i], branch.SubhaloMassType[i,4], branch.SubhaloSFR[i]
Output: 
135 45.3308 0.0147863
134 43.7982 0.00618325
133 43.8464 0.00186112
132 43.9853 0.000374569
131 44.4113 0.0
130 43.6553 0.0
129 43.3827 0.0125876
128 43.4569 0.00529478
127 43.8315 0.0531274
126 44.6463 0.0441715
...

32 0.00580271 0.921094
31 0.00488859 0.643712
30 0.00330759 0.392984
29 0.00213915 0.409696
28 0.00106858 0.427845
27 0.000718658 0.217374
26 7.60092e-05 0.0669792
25 0.0 0.0309659
24 0.0 0.0147244
23 0.0 0.00386036
22 0.0 0.00178407

Example 2

Print the total masses of the direct progenitors of a given subhalo. 
import readtreeHDF5
treedir = '/n/ghernquist/vrodrigu/MergerTrees/output/Subhalos/Illustris/L75n1820FP'
tree = readtreeHDF5.TreeDB(treedir)
snapnum = 85; subfind_id = 0
progs = tree.get_direct_progenitors(snapnum, subfind_id, keysel=['SubhaloMass'])
print progs.SubhaloMass
Output: 
[  5.46660938e+03   8.42722225e+00   5.96577488e-02   2.98375543e-02
   1.33459726e-02   1.75941736e-01   1.58780679e-01   2.00215373e-02
   1.11011732e-02   7.78871588e-03   1.76827796e-02   1.54313780e-02
   ...,
   8.81793071e-03   8.81793071e-03   8.81793071e-03   8.81793071e-03
   8.81793071e-03   8.81793071e-03   5.48702031e-02   2.28373185e-02
   1.89679563e-01   5.13658077e-02]

Example 3

For a given subhalo, count the number of progenitors at redshift z=3 (snapshot 60), each with total mass greater than 10^12 Msun/h. 
import readtreeHDF5
import numpy as np
treedir = '/n/ghernquist/vrodrigu/MergerTrees/output/Subhalos/Illustris/L75n1820FP'
tree = readtreeHDF5.TreeDB(treedir)
snapnum = 135; subfind_id = 0
subtree = tree.get_all_progenitors(snapnum, subfind_id, keysel=['SnapNum', 'SubhaloMass'])
bool_indices = np.logical_and(subtree.SnapNum == 60, subtree.SubhaloMass > 100.0)
print np.sum(bool_indices)
Output: 

5

Example 4

Print the gas fraction along the "future" branch of a subhalo, i.e. following the descendant links. 
import readtreeHDF5
import numpy as np
treedir = '/n/ghernquist/vrodrigu/MergerTrees/output/Subhalos/Illustris/L75n1820FP'
tree = readtreeHDF5.TreeDB(treedir)
snapnum = 60; subfind_id = 0
branch = tree.get_future_branch(snapnum, subfind_id, keysel=['SubhaloMassType'])
print branch.SubhaloMassType[:,0] / (branch.SubhaloMassType[:,0] + branch.SubhaloMassType[:,4])
Output: 
[ 0.81863177  0.817384    0.81635791  0.81488973  0.81344718  0.8120001
  0.80800343  0.80775136  0.81137735  0.80887872  0.81680244  0.81997812
  0.8095147   0.81191099  0.81993639  0.82624286  0.82041419  0.81541157
  ...,
  0.5013082   0.49303648  0.52275819  0.4496251   0.57909721  0.48233721
  0.58357298  0.65957606  0.68540359  0.71993554  0.72716904  0.77929348
  0.85557771  0.8851701   0.88132364]

Example 5

For a given subhalo, get all the subhalos from the same redshift that will eventually merge into the same object, and print their spatial location. 
import readtreeHDF5
import numpy as np
treedir = '/n/ghernquist/vrodrigu/MergerTrees/output/Subhalos/Illustris/L75n1820FP'
tree = readtreeHDF5.TreeDB(treedir)
snapnum = 60; subfind_id = 0
progs = tree.get_all_fellow_progenitors(snapnum, subfind_id, keysel=['SubhaloPos'])
print progs.SubhaloPos
Output: 
[[ 10846.50585938  50886.3203125   47641.7109375 ]
 [ 10823.27441406  50902.71484375  47633.87890625]
 [ 10791.67773438  50926.61328125  47626.51171875]
 ..., 
 [  9674.72363281  51196.58984375  47535.7421875 ]
 [ 11676.44042969  51291.84765625  48374.8984375 ]
 [  9742.10644531  50317.2109375   48392.4296875 ]]

Example 6

This last example does *not* make use of the reading script readtreeHDF5.py. Instead, it accesses the merger trees and offsets tables directly using PyTables. For a list of subhalos, we print the maximum DM mass they have ever had, along with the corresponding snapshot numbers. 
import numpy as np
import tables

# Constants
snapnum_last = 135
nsnaps = snapnum_last+1
parttype = 1  # DM particles

# List of subhalo indices at z=0
subfind_ids = np.arange(1000, dtype=np.int32)
num_subs = len(subfind_ids)

# Store maximum masses and snapshot numbers here:
max_masses = np.zeros(num_subs)
snap_numbers = np.zeros(num_subs, dtype=np.int32)

# Directory with the merger tree files
treedir = '/n/ghernquist/vrodrigu/MergerTrees/output/Subhalos/Illustris/L75n1820FP'

# Open merger tree and offsets files
f_tree = tables.openFile('%s/tree_extended.hdf5' % (treedir))
f_offsets = tables.openFile('%s/offsets/offsets_%s.hdf5' % (treedir, str(snapnum_last).zfill(3)))

# Get subhalo row numbers from offsets table,
# as well as the necessary info for extracting main branches
row_num = f_offsets.root.RowNum[:]
subhalo_id = f_offsets.root.SubhaloID[:]
main_leaf_progenitor_id = f_offsets.root.MainLeafProgenitorID[:]

# Now, for each subhalo in the list, get all the progenitors along
# the main branch and identify the most massive one.
for i in range(num_subs):
    # Get the masses and snapshot numbers along the main branch
    sub_index = subfind_ids[i]
    main_branch_length = main_leaf_progenitor_id[sub_index] - subhalo_id[sub_index] + 1
    locs = slice(row_num[sub_index], row_num[sub_index] + main_branch_length)
    masses = f_tree.root.SubhaloMassType[locs, parttype]
    snapnums = f_tree.root.SnapNum[locs]

    # Get max. mass and the corresponding snapshot number
    max_mass = 0.
    max_snapnum = -1
    for mass, snapnum in zip(masses, snapnums):
        if mass > max_mass:
            max_mass = mass
            max_snapnum = snapnum

    # Save to arrays
    max_masses[i] = max_mass
    snap_numbers[i] = max_snapnum

# Close files
f_tree.close()
f_offsets.close()

# Write results to file (masses in units of 10^10 Msun/h)
print '# subfind_id, max_mass [10^10 Msun/h], snapnum'
for i in range(num_subs):
    print '%d,%f,%d' % (subfind_ids[i], max_masses[i], snap_numbers[i])
Output: 
# subfind_id, max_mass (10^10 Msun/h), snapnum
0,20122.097656,135
1,2290.834961,131
2,2165.381348,119
3,336.076477,123
4,546.109863,100
5,296.441620,133
...
994,0.474405,118
995,10.420149,83
996,0.653409,102
997,0.378289,129
998,0.708962,106
999,0.441778,111

Example 7

Count the number of major mergers that the first 20 central subhalos have ever had (the script readsubfHDF5.py can be found at /n/home12/mvogelsberger/Python/OwnModules/). The mass ratio is determined by the maximum past stellar mass of the two progenitors. 
import numpy as np
import readtreeHDF5
import readsubfHDF5

def get_maximum_past_mass(subtree, rownum):
    # Get maximum past *stellar* mass of a subhalo
    parttype = 4
    branch_length = (subtree.MainLeafProgenitorID[rownum] -
                     subtree.SubhaloID[rownum] + 1)
    locs = slice(rownum, rownum + branch_length)
    masses = subtree.SubhaloMassType[locs, parttype]
    return np.max(masses)

def count_major_mergers(subtree, major_merger_ratio):
    num_major_mergers = 0
    root_sub_rownum = 0
    root_sub_id = subtree.SubhaloID[root_sub_rownum]

    # Iterate over "first progenitor" links (along main branch)
    first_prog_id = subtree.FirstProgenitorID[root_sub_rownum]
    while first_prog_id != -1:
        first_prog_rownum = root_sub_rownum + (first_prog_id - root_sub_id)
        first_prog_maxmass = get_maximum_past_mass(subtree, first_prog_rownum)

        # Iterate over "next progenitor" links
        next_prog_id = subtree.NextProgenitorID[first_prog_rownum]
        while next_prog_id != -1:
            next_prog_rownum = root_sub_rownum + (next_prog_id - root_sub_id)
            next_prog_maxmass = get_maximum_past_mass(subtree, next_prog_rownum)
            # Only if both masses are non-zero
            if first_prog_maxmass > 0 and next_prog_maxmass > 0:
                # Check if major meger (mass ratio determined by maximum past mass)
                if (next_prog_maxmass / first_prog_maxmass >= major_merger_ratio and
                    next_prog_maxmass / first_prog_maxmass <= 1.0/major_merger_ratio):
                    num_major_mergers += 1

            # Next iteration
            next_prog_id = subtree.NextProgenitorID[next_prog_rownum]

        # Next iteration
        first_prog_id = subtree.FirstProgenitorID[first_prog_rownum]

    return num_major_mergers

# -------------------------------- MAIN --------------------------------------

major_merger_ratio = 1.0/3.0
snapnum = 135
simulation = 'L75n1820FP'
tracking = 'Galaxies'  # Subhalos (DM) or Galaxies (baryons)
treedir = '/n/ghernquist/vrodrigu/MergerTrees/output/%s/Illustris/%s' % (tracking, simulation)
basedir = '/n/ghernquist/Illustris/Runs/%s/output' % (simulation)

# Look at the first few central subhalos
nsubs = 20
print 'Loading SUBFIND catalog...'
cat = readsubfHDF5.subfind_catalog(basedir, snapnum, long_ids=True, keysel=['GroupFirstSub'])
subfind_ids = [cat.GroupFirstSub[i] for i in range(nsubs)]

print 'Creating tree instance...'
tree = readtreeHDF5.TreeDB(treedir)

for i in range(nsubs):
    subtree = tree.get_all_progenitors(snapnum, subfind_ids[i],
            keysel=['SubhaloID', 'MainLeafProgenitorID', 'FirstProgenitorID',
            'NextProgenitorID', 'DescendantID', 'FirstSubhaloInFOFGroupID',
            'SubhaloMassType'])
    num_major_mergers = count_major_mergers(subtree, major_merger_ratio)
    print 'Subhalo %d had %d major mergers.' % (subfind_ids[i], num_major_mergers)
Output: 
Loading SUBFIND catalog...
Creating tree instance...
Subhalo 0 had 2 major mergers.
Subhalo 16937 had 2 major mergers.
Subhalo 30430 had 2 major mergers.
Subhalo 41088 had 2 major mergers.
Subhalo 51811 had 1 major mergers.
Subhalo 59384 had 1 major mergers.
Subhalo 66080 had 4 major mergers.
Subhalo 73663 had 6 major mergers.
Subhalo 80734 had 3 major mergers.
Subhalo 86186 had 2 major mergers.
Subhalo 93165 had 3 major mergers.
Subhalo 99148 had 4 major mergers.
Subhalo 104798 had 0 major mergers.
Subhalo 110566 had 2 major mergers.
Subhalo 114300 had 1 major mergers.
Subhalo 117343 had 2 major mergers.
Subhalo 120615 had 2 major mergers.
Subhalo 123773 had 1 major mergers.
Subhalo 127228 had 5 major mergers.
Subhalo 129770 had 2 major mergers.