Plotting number of halos per mass bin

Alexandres Lazar
  • 29 Apr '16

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.

Dylan Nelson
  • 29 Apr '16

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}$ ]')

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]
    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
Alexandres Lazar
  • 3 May '16

Awesome. Thank you.

  • Page 1 of 1