def get_halos_per_forest(f_in, Snap_Keys, haloID_field='ID', forestID_field
    ='ForestID'):
    """
    Determines the number of halos in each forest.

    The resulting Dictionary is nested with the outer-key given by the ForestID
    and the inner-key given by the snapshot field name.

    We also generate the offset for each Forest at each snapshot.  This is
    necessary because whilst Forest 5 may saved first at snapshot 20, it isn't
    necessarily saved first at snapshot 21.

    .. note::
        The default parameters are chosen to match the ASTRO3D Genesis trees as
        produced by VELOCIraptor + Treefrog.    

    Parameters
    ----------

    f_in: Open HDF5 file. 
        HDF5 file that contains the sorted trees.

    Snap_Keys: List of strings.
        List of keys that correspond to the fields containing the snapshot
        data.

    haloID_field: String. Default: 'ID'.
        Field name within the HDF5 file that corresponds to the unique halo ID.

    forestID_field: String. Default: 'ForestID'.
        Field name within the HDF5 file that corresponds to forest ID.
 
    Returns
    ----------

    NHalos_forest: Nested Dictionary.
        Nested dictionary that contains the number of halos for each Forest at
        each snapshot.  Outer-key is the ForestID and inner-key is the snapshot
        key.

    NHalos_forest_offset: Nested Dictionary.
        Nested dictionary that contains the offset for each Forest at each
        snapshot. Outer-key is the ForestID and inner-key is the snapshot key.
        This is required because whilst the tree is sorted by ForestID, the
        relative position of the tree can change from snapshot to snapshot.
    """
    import numpy as np
    import time
    import os
    from tqdm import tqdm
    start_time = time.time()
    print('')
    print(
        'Generating the dictionary for the number of halos in each forest at each snapshot.'
        )
    NHalos_forest = {}
    NHalos_forest_offset = {}
    for count, snap_key in enumerate(tqdm(Snap_Keys)):
        if len(f_in[snap_key][haloID_field]) == 0:
            continue
        halos_counted = 0
        halo_forestids = f_in[snap_key][forestID_field][:]
        forests_binned = np.bincount(halo_forestids)
        forestIDs = np.nonzero(forests_binned)[0]
        for forest_id in forestIDs:
            this_snap_NHalos = forests_binned[forest_id]
            try:
                NHalos_forest[forest_id][snap_key] = this_snap_NHalos
                NHalos_forest_offset[forest_id][snap_key] = halos_counted
            except KeyError:
                NHalos_forest[forest_id] = {snap_key: this_snap_NHalos}
                NHalos_forest_offset[forest_id] = {snap_key: halos_counted}
            halos_counted += this_snap_NHalos
    end_time = time.time()
    print('Creation of number of halos per forest took {0:3f} seconds.'.
        format(end_time - start_time))
    print('')
    return NHalos_forest, NHalos_forest_offset
