Public Data Access Overview / Discussion Forum

Hi,
I have been working on connecting particle information to the subhalo that they belong to, but I've run into a snag. I found a previous answer to a similar question with a similar function helpful, and while the code I've written worked previously for the snapshots I tested it on, it doesn't for some others in the same run. In particular, I am looking at black hole particles in TNG300-3, and while my code gave me good results for snapshot 99, using it for snapshot 11 has resulted in the data Subhalo/SnapByType from the offsets file outputting a different number of subhalos to SubhaloLenType from the group catalog.

I've been over my code several times, but I'm struggling to tell if this is a problem with the data I'm accessing here or my own code.

Many thanks.

Hi Kate,

Can you post a "minimal working example" of code which shows the problem?

For TNG300-3 at snapshot 11 there are 6142 black holes, of which 5440 are in groups (GroupLenType[:,5].sum()) and 5439 are in subhalos (SubhaloLenType[:,5].sum()).

Hi Dylan,

This is the code which is having the problem:

gName = 'Subhalo'
snap_n = 11

gcPath = basePath + '/groups_%03d/' % snap_n

with h5py.File(offsetPath(basePath,snap_n),'r') as f:
        offset_load = f['Subhalo/SnapByType']
        offsets = np.empty(np.shape(offset_load))
        offsets[:,:] = offset_load[:,:] 

ptType = 5
offsets_p = offsets[:,ptType]

for i in range(sim_metadata['num_files_groupcat']):
        with h5py.File(basePath+'output/groups_0{}/fof_subhalo_tab_0{}.{}.hdf5'.format(snap_n,snap_n,i),'r') as f:
                if i==0:
                        gcLenType = f['Subhalo/SubhaloLenType/'][:,ptType]
                else:
                        gcLenType = np.append(gcLenType,f['Subhalo/SubhaloLenType/'][:,ptType])

gcOffsetsType = offsets[:,ptType]

val = np.searchsorted(gcOffsetsType-1,p_index)
val = np.array([val - 1])
val = val.astype('int32')

flagFuzz = True
if flagFuzz:
        gcOffsetsMax = gcOffsetsType + gcLenType - 1
        ww = np.where(p_index > gcOffsetsMax[val])[0]

        if len(ww):
                val[ww] = -1

print('val: \n',val)

Hi Kate,

I suggest the following cleanups to the code, there isn't any need to do the loading so manually.

import illustris_python as il
import h5py
import numpy as np

def map_particle_index_to_parent_subhalo_index(p_index, flagFuzz=True):
    """ If flagFuzz is true, then the return has value -1 if the particle is 
    located in the FoF fuzz (outside any subhalo)."""

    basePath = 'sims.TNG/TNG300-3/output/'
    gName = 'Subhalo'
    snap_n = 11
    ptType = 5

    # load
    with h5py.File(il.groupcat.offsetPath(basePath,snap_n),'r') as f:
            gcOffsetsType = f['Subhalo/SnapByType'][:,ptType]

    gcLenType = il.groupcat.loadSubhalos(basePath, snap_n, 'SubhaloLenType')
    gcLenType = gcLenType[:,ptType]

    # find largest offset which is not past p_index
    val = np.searchsorted(gcOffsetsType-1,p_index)
    val = np.array([val - 1])
    val = val.astype('int32')

    if flagFuzz:
            gcOffsetsMax = gcOffsetsType + gcLenType - 1
            ww = np.where(p_index > gcOffsetsMax[val])[0]

            if len(ww):
                    val[ww] = -1

    return val

This runs fine for me and returns a reasonable number for the values input p_index values I tried. What input value produces an issue?