<p>Your approach is completely correct, just note three things:</p>

<h3>Logarithmic distribution</h3>

<p>First, since the distribution of masses is logarithmic in nature (as is most other things), be sure to bin them logarithmically. Otherwise you will oversample (undersample) the bins at the low-(high-)mass end.</p>

<h3>Comoving densities</h3>

<p>Second, to be able to compare mass functions at different redshifts, ones uses the <em>comoving</em> volume rather than the physical volume, such that the expansion of the Universe is factored out. The two are related as <span class="math-container">$V_\mathrm{com} = V_\mathrm{phys}(1+z)^3$</span>.</p>

<h3>Damn you, little <span class="math-container">$h$</span>!</h3>

<p>Finally, observers and modellers tend to use a slightly different definition of the unit volume. Whereas observers usually use <span class="math-container">$\mathrm{Mpc}$</span> for distances, and hence <span class="math-container">$\mathrm{Mpc}^{-3}$</span> for number densities, if your galaxies come from a cosmological simulation where the cosmological parameters can when tuned at will, it is custom to factor out the Hubble constant <span class="math-container">$H_0$</span>. In simulations, masses and distances are then measured in <span class="math-container">$h^{-1}M_\odot$</span> and <span class="math-container">$h^{-1}\mathrm{Mpc}$</span>, respectively, so number densities are measured in <span class="math-container">$h^3  \mathrm{Mpc}^{-3}$</span>. Here <span class="math-container">$h\equiv H_0\,/\,(100\,\mathrm{km}\,\mathrm{s}^{-1}\,\mathrm{Mpc}^{-1})$</span>.</p>

<p>This is probably a reminiscence of the time when the Hubble constant was rather uncertain. Nowadays, in my opinion there is no need to do this, but since everybody does, it's difficult to go against the stream. For a discussion of this issue, see <a href="https://arxiv.org/abs/1308.4150" rel="nofollow noreferrer">Croton (2013)</a>.</p>

<h3>Python code</h3>

<p>Since you've tagged the question with <a href="/questions/tagged/python" class="post-tag" title="show questions tagged &#39;python&#39;" rel="tag">python</a>, I wrote this little snippet that should do the work (I randomly chose <span class="math-container">$10^5\,\mathrm{Mpc}$</span> as your survey volume; note also that in this example I don't factor out <span class="math-container">$h$</span>):</p>

<pre><code>import numpy as np
import matplotlib.pyplot as plt

M     = np.loadtxt('Mstar.dat')         #Read stellar masses in Msun
logM  = np.log10(M)                     #Take logarithm
nbins = 10                              #Number of bins to divide data into
V     = 1e5                             #Survey volume in Mpc3
Phi,edg = np.histogram(logM,bins=nbins) #Unnormalized histogram and bin edges
dM    = edg[1] - edg[0]                 #Bin size
Max   = edg[0:-1] + dM/2.               #Mass axis
Phi   = Phi / V / dM                    #Normalize to volume and bin size

plt.clf()
plt.yscale('log')
plt.xlabel(r'<span class="math-container">$\log(M_\star\,/\,M_\odot)$</span>')
plt.ylabel(r'<span class="math-container">$\Phi\,/\,\mathrm{dex}^{-1}\,\mathrm{Mpc}^{-3}$</span>')
plt.plot(Max,Phi,ls='steps-post')
</code></pre>

CC BY-SA 3.0
2017-06-04T20:45:35.567Z
This is amazingly clear and helpful -- thank you very much! Especially for the python snippet!
quantumflash
CC BY-SA 3.0
2017-06-05T10:04:00.830Z
You're very welcome!
pela
CC BY-SA 4.0
2017-05-30T08:47:01.983Z
2019-03-26T19:24:30.640Z
2019-03-26T19:24:30.640Z
pela
<p>Suppose that I have an array of stellar masses of tens of thousands of galaxies, and that I also already know what the total "survey volume" is (in Mpc^3) within which all of those galaxies are contained. For context, my data is from a simulation of galaxy formation at redshift ~ 0, rather than from an observational simple -- so in principle I don't think I need to worry about incompleteness corrections at the faint (low-mass) end. </p>

<p>Given my array of stellar masses and the total survey volume in Mpc^3, how do I create and plot the "galaxy stellar mass function" that I always see people talking about? Isn't it just basically a histogram of the stellar mass array, and then divide the # of galaxies in every bin by the survey volume (it's almost just like a normalization constant)? </p>

<p>(Of course this question could equally apply to galaxy luminosity functions.) </p>

CC BY-SA 3.0
2017-05-29T10:03:54.137Z
Take a look at https://arxiv.org/pdf/1701.02734.pdf Chapt. 5
Py-ser
<galaxy><cosmology><mass><python><sky-survey>
How do I create a galaxy stellar mass function?
quantumflash
# How do I create a galaxy stellar mass function?

(Asked by: username_0 on 2017-05-29)

Suppose that I have an array of stellar masses of tens of thousands of galaxies, and that I also already know what the total "survey volume" is (in Mpc^3) within which all of those galaxies are contained. For context, my data is from a simulation of galaxy formation at redshift ~ 0, rather than from an observational simple -- so in principle I don't think I need to worry about incompleteness corrections at the faint (low-mass) end.

Given my array of stellar masses and the total survey volume in Mpc^3, how do I create and plot the "galaxy stellar mass function" that I always see people talking about? Isn't it just basically a histogram of the stellar mass array, and then divide the # of galaxies in every bin by the survey volume (it's almost just like a normalization constant)?

(Of course this question could equally apply to galaxy luminosity functions.)

**Tags:** galaxy, cosmology, mass, python, sky-survey

> **Comments:**
> **username_1**: Take a look at https://arxiv.org/pdf/1701.02734.pdf Chapt. 5  

---

## 1 answers

### Answer by username_2 on 2017-05-30. Score: 8 | Accepted answer

Your approach is completely correct, just note three things:
Logarithmic distribution
First, since the distribution of masses is logarithmic in nature (as is most other things), be sure to bin them logarithmically. Otherwise you will oversample (undersample) the bins at the low-(high-)mass end.
Comoving densities
Second, to be able to compare mass functions at different redshifts, ones uses the comoving volume rather than the physical volume, such that the expansion of the Universe is factored out. The two are related as $V_\mathrm{com} = V_\mathrm{phys}(1+z)^3$.
Damn you, little $h$!
Finally, observers and modellers tend to use a slightly different definition of the unit volume. Whereas observers usually use $\mathrm{Mpc}$ for distances, and hence $\mathrm{Mpc}^{-3}$ for number densities, if your galaxies come from a cosmological simulation where the cosmological parameters can when tuned at will, it is custom to factor out the Hubble constant $H_0$. In simulations, masses and distances are then measured in $h^{-1}M_\odot$ and $h^{-1}\mathrm{Mpc}$, respectively, so number densities are measured in $h^3  \mathrm{Mpc}^{-3}$. Here $h\equiv H_0\,/\,(100\,\mathrm{km}\,\mathrm{s}^{-1}\,\mathrm{Mpc}^{-1})$.
This is probably a reminiscence of the time when the Hubble constant was rather uncertain. Nowadays, in my opinion there is no need to do this, but since everybody does, it's difficult to go against the stream. For a discussion of this issue, see Croton (2013).
Python code
Since you've tagged the question with python, I wrote this little snippet that should do the work (I randomly chose $10^5\,\mathrm{Mpc}$ as your survey volume; note also that in this example I don't factor out $h$):
import numpy as np
import matplotlib.pyplot as plt

M     = np.loadtxt('Mstar.dat')         #Read stellar masses in Msun
logM  = np.log10(M)                     #Take logarithm
nbins = 10                              #Number of bins to divide data into
V     = 1e5                             #Survey volume in Mpc3
Phi,edg = np.histogram(logM,bins=nbins) #Unnormalized histogram and bin edges
dM    = edg[1] - edg[0]                 #Bin size
Max   = edg[0:-1] + dM/2.               #Mass axis
Phi   = Phi / V / dM                    #Normalize to volume and bin size

plt.clf()
plt.yscale('log')
plt.xlabel(r'$\log(M_\star\,/\,M_\odot)$')
plt.ylabel(r'$\Phi\,/\,\mathrm{dex}^{-1}\,\mathrm{Mpc}^{-3}$')
plt.plot(Max,Phi,ls='steps-post')

> **Comments:**
> **username_0**: This is amazingly clear and helpful -- thank you very much! Especially for the python snippet!  
> **username_2**: You're very welcome!  

---

