Public Data Access Overview / Discussion Forum

Without using and any type of theoretical formalism, I'm wanting to plot the number of dark matter halos per logarithmic bin of mass, reminiscent that of the halo mass function.

I of course use the catalog of the number of dark matter halos, but how can i go about plotting the number of halos between mass M and M+dM? Would I use the size of the box as the volume and the number of mass to calculate the number density? I'm a bit lost how to approach and go about out this.

Hi Alexandres,

Below some snippets of code for a stellar mass function, should be very similar for a DM halo mass function.

ax.set_xlabel('Galaxy Stellar Mass [ log M$_{\\rm sun}$ ] [ < 1r$_{1/2}$, < 2r$_{1/2}$, or total ]')
ax.set_ylabel('$\Phi$ [ Mpc$^{-3}$ dex$^{-1}$ ]')
ax.set_yscale('log')

mts = ['SubhaloMassInRadType','SubhaloMassInHalfRadType','SubhaloMassType']
gc = groupCat(sP, fieldsHalos=['GroupFirstSub','Group_M_Crit200'], fieldsSubhalos=mts)

# centrals only?
if centralsOnly:
    wHalo = np.where((gc['halos']['GroupFirstSub'] >= 0))
    w = gc['halos']['GroupFirstSub'][wHalo]
else:
    w = np.arange(gc['subhalos']['count'])

for mt in mts:
    xx = gc['subhalos'][mt][w,partTypeNum('stars')]
    xx = codeMassToLogMsun(xx)
    xm, ym = running_histogram(xx, binSize=binSize, normFac=boxSizeCubicMpc*binSize)
    ax.plot(xm, ym)

At the high mass end (with low numbers) the binning decisions can become important, this is just one approach:

def running_histogram(X, nBins=100, binSize=None, normFac=None):
    """ Create a adaptive histogram of a (x) point set using some number of bins. """
    if binSize is not None:
        nBins = round( (X.max()-X.min()) / binSize )

    bins = np.linspace(X.min(),X.max(), nBins)
    delta = bins[1]-bins[0]

    running_h   = []
    bin_centers = []

    for i, bin in enumerate(bins):
        binMax = bin + delta
        w = np.where((X >= bin) & (X < binMax))

        if len(w[0]):
            running_h.append( len(w[0]) )
            bin_centers.append( np.nanmedian(X[w]) )

    if normFac is not None:
        running_h /= normFac

    return bin_centers, running_h

Awesome. Thank you.

Hi Dylan,

Hope all is well. I realized that the API signature is now different - this snippet was from 2016 - any chance you could past this snippet using the new ill.groupcat.loadHalos(...) syntax please?

Hi Andre,

The API and data are all unchanged since forever. The code above is just to give an idea, it won't run as is. You can adapt it into what you need, i.e. just replacing the gc = groupCat() line with a halos = il.groupcat.loadHalos() and subhalos = il.groupcat.loadSubhalos().