diff --git a/skypy/galaxies/redshift.py b/skypy/galaxies/redshift.py
index f048a357..ecf64d33 100644
--- a/skypy/galaxies/redshift.py
+++ b/skypy/galaxies/redshift.py
@@ -126,7 +126,7 @@ def f(lnx, a):
         gam[i], _ = scipy.integrate.quad(f, lnxmin[i], np.inf, args=(alpha[i],))
 
     # comoving number density is normalisation times upper incomplete gamma
-    density = phi_star*gam
+    density = phi_star*gam*np.power(cosmology.h, 3)
 
     # sample redshifts from the comoving density
     return redshifts_from_comoving_density(redshift=redshift, density=density,
diff --git a/skypy/galaxies/tests/test_redshift.py b/skypy/galaxies/tests/test_redshift.py
index d41e92ac..a5208f9d 100644
--- a/skypy/galaxies/tests/test_redshift.py
+++ b/skypy/galaxies/tests/test_redshift.py
@@ -32,7 +32,7 @@ def test_schechter_lf_redshift():
     x_min = 10.**(-0.4*(M_lim - M_star))
 
     # density with factor from upper incomplete gamma function
-    density = phi_star*gamma(alpha+1)*gammaincc(alpha+1, x_min)
+    density = phi_star*gamma(alpha+1)*gammaincc(alpha+1, x_min)*np.power(cosmo.h, 3)
 
     # turn into galaxies/surface area
     density *= (sky_area * cosmo.differential_comoving_volume(z)).to_value('Mpc3')