Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Delta sigma from 2d #696

Merged
merged 6 commits into from Mar 9, 2017
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view

Large diffs are not rendered by default.

Expand Up @@ -6,32 +6,32 @@ Galaxy Catalog Analysis Example: Galaxy-galaxy lensing
In this example, we'll show how to calculate :math:`\Delta\Sigma(r),`
the galaxy-galaxy lensing signal of a mock catalog.

There is also an IPython Notebook in the following location that can be
There is also a Jupyter Notebook in the following location that can be
used as a companion to the material in this section of the tutorial:


**halotools/docs/notebooks/galcat_analysis/basic_examples/galaxy_catalog_analysis_tutorial3.ipynb**

By following this tutorial together with this notebook,
you can play around with your own variations of the calculation
as you learn the basic syntax.
By following this tutorial together with this notebook,
you can play around with your own variations of the calculation
as you learn the basic syntax.

Generate a mock galaxy catalog
Generate a mock galaxy catalog
---------------------------------
Let's start out by generating a mock galaxy catalog into an N-body
simulation in the usual way. Here we'll assume you have the *z=0*
rockstar halos for the bolshoi simulation, as this is the
default halo catalog.
default halo catalog.

.. code:: python

from halotools.empirical_models import PrebuiltSubhaloModelFactory
model = PrebuiltSubhaloModelFactory('smhm_binary_sfr')
model = PrebuiltSubhaloModelFactory('behroozi10')
from halotools.sim_manager import CachedHaloCatalog
halocat = CachedHaloCatalog(simname = 'bolshoi', redshift = 0, halo_finder = 'rockstar')
model.populate_mock(halocat)

Extract subsamples of galaxies and dark matter particles
Extract subsamples of galaxies and dark matter particles
------------------------------------------------------------------
Predictions for galaxy-galaxy lensing are calculated from the
cross-correlation between the galaxy positions and the dark matter
Expand All @@ -44,13 +44,13 @@ positions stored in the ``ptcl_table`` attribute of the mock.
py = model.mock.ptcl_table['y']
pz = model.mock.ptcl_table['z']

As described in :ref:`mock_obs_pos_formatting`,
functions in the `~halotools.mock_observables` package
such as `~halotools.mock_observables.delta_sigma` take array inputs in a
specific form: a (*Npts, 3)*-shape Numpy array. You can use the
`~halotools.mock_observables.return_xyz_formatted_array` convenience
function for this purpose, which has a built-in *mask* feature
that we'll also demonstrate to select a random downsampling of :math:`10^{5}`
As described in :ref:`mock_obs_pos_formatting`,
functions in the `~halotools.mock_observables` package
such as `~halotools.mock_observables.delta_sigma` take array inputs in a
specific form: a (*Npts, 3)*-shape Numpy array. You can use the
`~halotools.mock_observables.return_xyz_formatted_array` convenience
function for this purpose, which has a built-in *mask* feature
that we'll also demonstrate to select a random downsampling of :math:`10^{5}`
dark matter particles.

.. code:: python
Expand All @@ -68,53 +68,66 @@ dark matter particles.
ptcl_mask = np.where(sorted_randoms < sorted_randoms[Nptcls_to_keep])[0]
particle_positions = return_xyz_formatted_array(px, py, pz, mask = ptcl_mask)

Now we'll extract the *x, y, z* positions of various subsamples of our galaxies.
Now we'll extract the *x, y, z* positions of various subsamples of our galaxies.

.. code:: python

x = model.mock.galaxy_table['x']
y = model.mock.galaxy_table['y']
z = model.mock.galaxy_table['z']

mstar11_mask = model.mock.galaxy_table['stellar_mass'] > 1e11
mstar11_positions = return_xyz_formatted_array(x, y, z, mask = mstar11_mask)

mstar105_mask = (model.mock.galaxy_table['stellar_mass'] > 10**10.3) & (model.mock.galaxy_table['stellar_mass'] < 10**10.7)
mstar105_positions = return_xyz_formatted_array(x, y, z, mask = mstar105_mask)

mstar105_central_mask = mstar105_mask * (model.mock.galaxy_table['halo_upid'] == -1)
mstar105_central_positions = return_xyz_formatted_array(x, y, z, mask = mstar105_central_mask)

mstar105_satellite_mask = mstar105_mask * (model.mock.galaxy_table['halo_upid'] != -1)
mstar105_satellite_positions = return_xyz_formatted_array(x, y, z, mask = mstar105_satellite_mask)

Calculate :math:`\Delta\Sigma(R_{\rm p})`
-------------------------------------------------------------

As of Halotools 0.5, the `~halotools.mock_observables.delta_sigma` function requires you to
specify the particle masses of your simulation, as well as the factory by which you have
randomly downsampled the particles in the snapshot to perform your calculation.

.. code:: python

from halotools.mock_observables import delta_sigma


rp_bins = np.logspace(-1, 1, 15)

particle_masses = halocat.particle_mass
period=model.mock.Lbox
downsampling_factor = (halocat.num_ptcl_per_dim**3)/float(len(particle_positions))

rp_bins = np.logspace(-1,1,15)
pi_max = 40

result_mstar11_in_mpc = delta_sigma(mstar11_positions, particle_positions,
rp_bins, pi_max, period=model.mock.Lbox)

result_mstar105_in_mpc = delta_sigma(mstar105_positions, particle_positions,
rp_bins, pi_max, period=model.mock.Lbox)

result_mstar105_central_in_mpc = delta_sigma(mstar105_central_positions, particle_positions,
rp_bins, pi_max, period=model.mock.Lbox)

result_mstar105_satellite_in_mpc = delta_sigma(mstar105_satellite_positions, particle_positions,
rp_bins, pi_max, period=model.mock.Lbox)


Recall that all Halotools length units are comoving and in Mpc/h. However, the conventional units to
plot :math:`\Delta\Sigma` are :math:`h*M_{\odot}/pc^2`, since in those units the galaxy-galaxy
lensing signal is roughly order unity for typical :math:`L_{\ast}` galaxy samples.
So now we convert units and plot the results.

rp, result_mstar11_in_mpc = delta_sigma(mstar11_positions, particle_positions,
particle_masses, downsampling_factor,
rp_bins, period, cosmology=halocat.cosmology, num_threads='max')

rp, result_mstar105_in_mpc = delta_sigma(mstar105_positions, particle_positions,
particle_masses, downsampling_factor,
rp_bins, period, cosmology=halocat.cosmology, num_threads='max')

rp, result_mstar105_central_in_mpc = delta_sigma(mstar105_central_positions, particle_positions,
particle_masses, downsampling_factor,
rp_bins, period, cosmology=halocat.cosmology, num_threads='max')

rp, result_mstar105_satellite_in_mpc = delta_sigma(mstar105_satellite_positions, particle_positions,
particle_masses, downsampling_factor,
rp_bins, period, cosmology=halocat.cosmology, num_threads='max')


Recall that all Halotools length units are comoving and in Mpc/h. However, the conventional units to
plot :math:`\Delta\Sigma` are :math:`h*M_{\odot}/pc^2`, since in those units the galaxy-galaxy
lensing signal is roughly order unity for typical :math:`L_{\ast}` galaxy samples.
So now we convert units and plot the results.

.. code:: python

Expand All @@ -126,24 +139,27 @@ So now we convert units and plot the results.



Plot the results
Plot the results
~~~~~~~~~~~~~~~~~~~~
.. code:: python

plt.plot(rp_bins, result_mstar11_in_pc, label=r'All galaxies: $M_{\ast} > 10^{11}M_{\odot}$')
plt.plot(rp_bins, result_mstar105_in_pc, label=r'All galaxies: $M_{\ast} \approx 10^{10.5}M_{\odot}$')
plt.plot(rp_bins, result_mstar105_satellite_in_pc, label=r'Satellites: $M_{\ast} \approx 10^{10.5}M_{\odot}$')
plt.plot(rp_bins, result_mstar105_central_in_pc, label=r'Centrals: $M_{\ast} \approx 10^{10.5}M_{\odot}$')

plt.xlim(xmin = 0.1, xmax = 10)
plt.ylim(ymin = 0.1, ymax = 50000)
plt.loglog()
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.xlabel(r'$R_{\rm p} $ $\rm{[Mpc / h]}$', fontsize=22)
plt.ylabel(r'$\Delta\Sigma(R_{\rm p})$ $[h M_{\odot} / {\rm pc}^2]$', fontsize=22)
plt.legend(loc='best', fontsize=16)
fig, ax = plt.subplots(1, 1, figsize=(8, 6))

__=plt.loglog()

__=ax.plot(rp, result_mstar11_in_pc, label=r'All galaxies: $M_{\ast} > 10^{11}M_{\odot}$')
__=ax.plot(rp, result_mstar105_satellite_in_pc, label=r'Satellites: $M_{\ast} \approx 10^{10.5}M_{\odot}$')
__=ax.plot(rp, result_mstar105_in_pc, label=r'All galaxies: $M_{\ast} \approx 10^{10.5}M_{\odot}$')
__=ax.plot(rp, result_mstar105_central_in_pc, label=r'Centrals: $M_{\ast} \approx 10^{10.5}M_{\odot}$')

__=ax.set_xlim(xmin = 0.1, xmax = 10)
__=ax.set_ylim(ymin = 0.5, ymax = 200)

__=ax.set_xlabel(r'$R_{\rm p} $ $\rm{[Mpc / h]}$', fontsize=16)
__=ax.set_ylabel(r'$\Delta\Sigma(R_{\rm p})$ $[h M_{\odot} / {\rm pc}^2]$', fontsize=16)
__=ax.legend(loc='best', fontsize=13)
__=plt.xticks(fontsize=15); plt.yticks(fontsize=15)

.. image:: gg_lensing_tutorial3.png

This tutorial continues with :ref:`galaxy_catalog_analysis_tutorial4`.
This tutorial continues with :ref:`galaxy_catalog_analysis_tutorial4`.
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
1 change: 1 addition & 0 deletions halotools/mock_observables/__init__.py
Expand Up @@ -18,3 +18,4 @@
from .large_scale_density import *
from .counts_in_cells import *
from .occupation_stats import hod_from_mock, get_haloprop_of_galaxies
from .surface_density import *
1 change: 1 addition & 0 deletions halotools/mock_observables/surface_density/__init__.py
@@ -0,0 +1 @@
from .delta_sigma import delta_sigma