Author: Julien Peloton @JulienPeloton
Last Verifed to Run: 2019-01-04
Estimated running time: < 5 min.
This notebook, with the help of Apache Spark, inspects the peaks found in the halo mass distribution in the companion notebook.
Summary:
If we look at the distribution of halo masses in cosmoDC2, all seem OK. But if we look at the same distribution after filtering objects according to the stellar mass, we start to see peaks in the distribution. These peaks are regularly spaced in log, and their position is independent on the stellar mass cut applied.
Useful links:
Logistics:
This notebook is intended to be run through the JupyterHub NERSC interface with the desc-pyspark kernel. The kernel is automatically installed in your environment when you use the kernel setup script:
source /global/common/software/lsst/common/miniconda/kernels/setup.sh
import numpy as np
import matplotlib
%matplotlib inline
import matplotlib.pyplot as plt
# Path to the data
fn = "/global/cscratch1/sd/peloton/cosmodc2/xyz_v1.1.4_mass_and_mag.parquet"
# Load the data - this a lazy operation, no data movement yet!
df = spark.read.format("parquet").load(fn)
# Let's inspect the schema
df.printSchema()
# Number of objects in the catalog
print("Number of rows: {}".format(df.count()))
# We reject synthetic halos - see for example discussion in
# https://github.com/LSSTDESC/cosmodc2/issues/82
# In addition we select halo members according to their stellar mass
stellar_masses_cut = [1e9, 1e10, 5e10, 1e11]
full_data = []
for stellar_mass_cut in stellar_masses_cut:
df_filt = df.filter("halo_id > 0").filter("stellar_mass > {}".format(stellar_mass_cut))
# Group data by halos and compute the number of halo members
df_disp = df_filt.groupBy("halo_id").count()
# Add back the original DataFrame columns
# and select only central member for halo
# (unique halo mass for a halo)
data_joined = df_filt.join(df_disp, "halo_id")\
.filter("is_central == True")\
.select("halo_mass", 'count')\
.dropna()
# Collect the data from the executors to the driver
data = data_joined.collect()
full_data.append(data)
(for the version without cut, see the main notebook)
matplotlib.rcParams.update({'font.size': 17})
fig = plt.figure(figsize=(10, 7))
plt.title("Histogram of halo masses ($n_{{gal}}/halo > 0$)")
# Plot peak locations (empirical!)
for k in range(90):
plt.axvline(k*0.15, ls='--', color='k', alpha=0.5)
# Plot halo mass data
for index, stellar_mass_cut in enumerate(stellar_masses_cut):
mass, count = np.transpose(full_data[index])
print(r"Number of entries ($M_*$ > {:.1e} $M_\odot$): {}".format(stellar_mass_cut, len(mass)))
plt.hist(np.log10(mass), bins=200, alpha=0.5, label='$M_*$ > {:.1e} $M_\odot$'.format(stellar_mass_cut))
plt.xlim(10.3, 15)
plt.xlabel(r'$\log({\rm M}_h \, [{\rm M}_\odot])$')
plt.yscale('log')
plt.legend(title="Only central galaxies with:")
plt.show()
The position $p$ of the peaks is independent of the cut on the stellar mass, and it seems to follow:
$$ \log(M_{h}^{p}) = 0.15 * p $$Let's plot the same distribution, but selecting only halos with at least 2 members
matplotlib.rcParams.update({'font.size': 17})
mincount = 1
fig = plt.figure(figsize=(10, 7))
plt.title("Histogram of halo masses ($n_{{gal}}/halo >$ {})".format(mincount))
# Plot peak locations (empirical!)
for k in range(90):
plt.axvline(k*0.15, ls='--', color='k', alpha=0.5)
# Plot halo mass data
for index, stellar_mass_cut in enumerate(stellar_masses_cut):
mass, count = np.transpose(full_data[index])
mask = count > mincount
print(r"Number of entries ($M_*$ > {:.1e} $M_\odot$): {}".format(stellar_mass_cut, len(mass[mask])))
plt.hist(np.log10(mass[mask]), bins=200, alpha=0.5, label='$M_*$ > {:.1e} $M_\odot$'.format(stellar_mass_cut))
plt.xlim(10.3, 15)
plt.xlabel(r'$\log({\rm M}_h \, [{\rm M}_\odot])$')
plt.yscale('log')
plt.legend(title="Only central galaxies with:")
plt.show()
The peaks are even sharper at low mass.
Question: Is that expected?
Answer: The saw tooth shape is indeed expected (but probably not wanted!) from the way N-body simulations are done. A detailed discussion can be found at https://github.com/LSSTDESC/DC2-analysis/issues/57, and an action item for the simulation has been opened at https://github.com/LSSTDESC/cosmodc2/issues/84.
It has been suggested to not use the stellar mass for selecting halos (not a direct observable), but rather a more physically motivated quantity like magnitudes/colors. In this part, we select the apparent magnitude, lensed, in i band. This might not be the best quantity to look at, but just want to play with something else than the stellar mass.
# let's select halo masses, and magnitudes (i, g, u) for
# the central galaxy.
data = df.filter("halo_id > 0")\
.filter("is_central == True")\
.select(["halo_mass", "mag_i", "mag_g", "mag_u"])\
.sample(fraction=0.05).collect()
mass, mag_i, mag_g, mag_u = np.transpose(data)
print("Number of rows: {}".format(len(mass)))
Let's have a look at the 2D histogram halo mass - magnitude in i band:
import seaborn as sns
joint_kws = dict(gridsize=200, mincnt=1)
g = sns.jointplot(
np.log10(mass),
mag_i,
height=8, space=0.5,
kind='hex', color='k',
xlim=(10.5, 12), ylim=(20, 35),
joint_kws=joint_kws,
marginal_kws=dict(bins=200, rug=True))
g.set_axis_labels(
r'$\log({\rm M}_h \, [{\rm M}_\odot])$',
'Apparent magnitude, lensed, in i band (lsst)')
plt.show()
Hum, the stripes are still there... Let's have a look at 1D halo mass distribution as a function of magnitude:
# We reject synthetic halos - see for example discussion in
# https://github.com/LSSTDESC/cosmodc2/issues/82
# In addition we select halo members according to their stellar mass
stellar_masses_cut = [25.0, 22.5, 20.0]
full_data = []
for stellar_mass_cut in stellar_masses_cut:
df_filt = df.filter("halo_id > 0").filter("mag_i < {}".format(stellar_mass_cut))
# Group data by halos and compute the number of halo members
df_disp = df_filt.groupBy("halo_id").count()
# Add back the original DataFrame columns
# and select only central member for halo
# (unique halo mass for a halo)
data_joined = df_filt.join(df_disp, "halo_id")\
.filter("is_central == True")\
.select("halo_mass", 'count')\
.dropna()
# Collect the data from the executors to the driver
data = data_joined.collect()
full_data.append(data)
matplotlib.rcParams.update({'font.size': 17})
fig = plt.figure(figsize=(10, 7))
plt.title("Histogram of halo masses ($n_{{gal}}/halo > 0$)")
# Plot peak locations (empirical!)
for k in range(90):
plt.axvline(k*0.15, ls='--', color='k', alpha=0.5)
# Plot halo mass data
for index, stellar_mass_cut in enumerate(stellar_masses_cut):
mass, count = np.transpose(full_data[index])
print(r"Number of entries (mag_i_lsst < {}): {}".format(stellar_mass_cut, len(mass)))
plt.hist(np.log10(mass), bins=200, alpha=0.5, label='mag_i_lsst < {}'.format(stellar_mass_cut))
plt.xlim(10.3, 15)
plt.xlabel(r'$\log({\rm M}_h \, [{\rm M}_\odot])$')
plt.yscale('log')
plt.legend(title="Only central galaxies with:")
plt.show()
Similar to what is seen for cut on the stellar mass... According to Andrew Hearin in #57, this makes sense as
... in the cosmoDC2 model, restframe flux derives from stellar mass,
which in turn drives from halo mass, so it is expected that this discreteness
propagates through to these other variables.