Example Applications

Below are some examples demonstrating interesting combinations of ytree functionality. Each of the scripts shown below can be found in the doc/source/examples directory. If you have made something not seen here, please considering adding it to this document.

Plot the Tree of the Most Massive Halo

Script: plot_most_massive.py

Below we make a plot of the most massive halo in the arbor. We use the NumPy argmax function to get the index within the arbor of the most massive halo.

import ytree

a = ytree.load("consistent_trees/tree_0_0_0.dat")

imax = a["mass"].argmax()
my_tree = a[imax]
print(f"Most massive halo is {my_tree} with M = {my_tree['mass']}.")

p = ytree.TreePlot(my_tree)
p.min_mass_ratio = 0.001
p.save("most_massive.png")

We use the min_mass_ratio attribute to plot only halos with masses of at least 10^-3 of the main halo.

Plot the Tree with the Most Halos

Script: plot_most_halos.py

Similar to above, it is often useful to find the tree containing the most halos. To do this, we make an array containing the sizes of all trees using the tree_size attribute of the TreeNode class. The Arbor class’s arr method is useful for creating unyt_array objects with the unit system of the dataset.

import ytree

a = ytree.load("consistent_trees/tree_0_0_0.dat")

tree_size = a.arr([t.tree_size for t in a])
imax = tree_size.argmax()
my_tree = a[imax]
print(f"Tree with most halos is {my_tree} with {my_tree.tree_size} halos.")

p = ytree.TreePlot(my_tree)
p.min_mass_ratio = 0.001
p.save("most_halos.png")

Halo Age (a50)

Script: halo_age.py

Note

This script includes extra code to make it run within the test suite. To run conventionally, remove the lines indicated in the header of script.

One way to define the age of a halo is by calculating the scale factor when it reached 50% of its current mass. This is often referred to as “a50”. In the example below, this is calculated by linearly interpolating from the mass of the main progenitor.

def calc_a50(node):
    # main progenitor masses
    pmass = node["prog", "mass"]

    mh = 0.5 * node["mass"]
    m50 = pmass <= mh

    if not m50.any():
        th = node["scale_factor"]
    else:
        pscale = node["prog", "scale_factor"]
        # linearly interpolate
        i = np.where(m50)[0][0]
        slope = (pscale[i - 1] - pscale[i]) / (pmass[i - 1] - pmass[i])
        th = slope * (mh - pmass[i]) + pscale[i]

    node["a50"] = th

Then, we setup an Analysis Pipeline including this function and use parallel_nodes to loop over all halos in the dataset in parallel.

        a = ytree.load("tiny_ctrees/locations.dat")
        a.add_analysis_field("a50", "")

        ap = ytree.AnalysisPipeline()
        ap.add_operation(calc_a50)

        trees = list(a[:])
        for tree in ytree.parallel_nodes(trees, filename="halo_age"):
            yt.mylog.info(f"Processing {tree}.")
            ap.process_target(tree)

Finally, we reload the saved data and print the age of the first halo.

        if yt.is_root():
            a2 = ytree.load("halo_age/halo_age.h5")
            print(a2[0]["a50"])

Do the following to run the script on two processors:

$ mpirun -np 2 python halo_age.py

Significance

Script: halo_significance.py

Note

This script includes extra code to make it run within the test suite. To run conventionally, remove the lines indicated in the header of script.

Brought to you by John Wise, a halo’s significance is calculated by recursively summing over all ancestors the mass multiplied by the time between snapshots. When determining the main progenitor of a halo, the significance measure will select for the ancestor with the deeper history instead of just the higher mass. This can be helpful in cases of near 1:1 mergers.

First, we define a function that calculates the significance for every halo in a single tree.

def calc_significance(node):
    if node.descendent is None:
        dt = 0.0 * node["time"]
    else:
        dt = node.descendent["time"] - node["time"]

    sig = node["mass"] * dt
    if node.ancestors is not None:
        for anc in node.ancestors:
            sig += calc_significance(anc)

    node["significance"] = sig
    return sig

Then, we use the The AnalysisPipeline to calculate the significance for all trees and save a new dataset. Because the calc_significance function defined above works on all halos in a given tree at once, we parallelize this by allocating a whole tree to each processor using the parallel_trees function.

        a = ytree.load("tiny_ctrees/locations.dat")
        a.add_analysis_field("significance", "Msun*Myr")

        ap = ytree.AnalysisPipeline()
        ap.add_operation(calc_significance)

        trees = list(a[:])
        for tree in ytree.parallel_trees(trees, filename="halo_significance"):
            yt.mylog.info(f"Processing {tree}.")
            ap.process_target(tree)

After loading the new arbor, we use the set_selector function to use the new significance field to determine the progenitor line.

        if yt.is_root():
            a2 = ytree.load("halo_significance/halo_significance.h5")
            a2.set_selector("max_field_value", "significance")
            prog = list(a2[0]["prog"])
            print(prog)

Do the following to run the script on two processors:

$ mpirun -np 2 python halo_significance.py

Merger and Smooth Accretion Rates

Script: merger_accretion_rates.py

Note

This script includes extra code to make it run within the test suite. To run conventionally, remove the lines indicated in the header of script.

The growth of a dark matter halo can be decomposed into two components: mass gained through mergers with other halos, and mass gained through smooth accretion of diffuse material. Below we define a function that calculates both rates for each halo node.

The merger rate is calculated by summing the masses of all direct ancestors, subtracting the main progenitor mass, and dividing by the time between them. The smooth accretion rate is then inferred as the residual between the total mass growth rate and the merger rate.

def calc_accretion_rates(node):
    ancestors = list(node.ancestors)
    if not ancestors:
        return

    main_prog = ancestors[np.argmax([p["mass"] for p in ancestors])]

    dt = node["time"] - main_prog["time"]
    if dt == 0:
        return

    dm_total = node["mass"] - main_prog["mass"]
    m_all_progs = sum(p["mass"] for p in ancestors)
    dm_merger = m_all_progs - main_prog["mass"]

    node["mdot_merger"] = dm_merger / dt
    node["mdot_accretion"] = (dm_total - dm_merger) / dt

Then, we setup an Analysis Pipeline and use parallel_nodes to loop over all halos in parallel.

        a = ytree.load("tiny_ctrees/locations.dat")
        a.add_analysis_field("mdot_merger", units="Msun/yr", default=0.0)
        a.add_analysis_field("mdot_accretion", units="Msun/yr", default=0.0)

        ap = ytree.AnalysisPipeline()
        ap.add_operation(calc_accretion_rates)

        trees = list(a[:])
        for tree in ytree.parallel_nodes(trees, filename="accretion_rates"):
            yt.mylog.info(f"Processing {tree}.")
            ap.process_target(tree)

Do the following to run the script on two processors:

$ mpirun -np 2 python merger_accretion_rates.py