Generate mock data for a cluster ensemble

Generate cluster ensemble with random masses and redshifts

%load_ext autoreload
%autoreload 2

import sys
import os
import numpy as np
from astropy.table import Table
from numpy import random
import scipy
import matplotlib.pyplot as plt
import clmm
from clmm import GalaxyCluster, ClusterEnsemble, GCData
from clmm import Cosmology
from clmm.support import mock_data as mock

clmm.__version__

'1.12.0'

np.random.seed(11)

cosmo = Cosmology(H0=71.0, Omega_dm0=0.265 - 0.0448, Omega_b0=0.0448, Omega_k0=0.0)

Generating a cluster catalog and associated source catalogs

Below, we randomly generate the masses, redshifts, concentrations and coordinates of an ensemble of n_clusters clusters. For simplicity, the drawing is performed uniformly in logm and redshift (instead of following a halo mass function).

# redshift and mass range of the galaxy clusters
z_bin = [0.2, 0.25]
logm_bin = np.array([14, 14.1])  # Solar Mass

# number of clusters in the ensemble
n_clusters = 30

# random draw in the mass and redshift range (for simplicity, uniform instead of following an actual mass function)
cluster_m = 10 ** (
    (logm_bin[1] - logm_bin[0]) * np.random.random(n_clusters) + logm_bin[0]
)  # in M_sun
cluster_z = (z_bin[1] - z_bin[0]) * np.random.random(n_clusters) + z_bin[0]

# random normal draw of cluster concentration, around c_mean
c_mean = 4.0
lnc = abs(np.log(c_mean) + 0.01 * np.random.randn(n_clusters))
concentration = np.exp(lnc)

# randomly draw cluster positions on the sky
ra = np.random.random(n_clusters) * 360  # from 0 to 360 deg
sindec = np.random.random(n_clusters) * 2 - 1
dec = np.arcsin(sindec) * 180 / np.pi  # from -90 to 90 deg

Background galaxy catalog generation

For each cluster of the ensemble, we use mock_data to generate a background galaxy catalog and store the results in a GalaxyCluster object. Note that: - The cluster density profiles follow the NFW parametrisation; - The source redshifts follow the Chang et al. distribution and have associated pdfs; - The shapes include shape noise and shape measurement errors; - Background galaxy catalogs are independent, even if the clusters are close (i.e., no common galaxy between two catalogs). - For each cluster we then compute - the tangential and cross \(\Delta\Sigma\) for each background galaxy - the weights w_ls to be used to compute the corresponding radial profiles (see demo_compute_deltasigma_weights.ipynb notebook for details)

The cluster objects are then stored in gclist.

import warnings

warnings.simplefilter(
    "ignore"
)  # just to prevent warning print out when looping over the cluster ensemble below

gclist = []
# number of galaxies in each cluster field (alternatively, can use the galaxy density instead)
n_gals = 10000
# ngal_density = 10

for i in range(n_clusters):
    # generate background galaxy catalog for cluster i
    noisy_data_z = mock.generate_galaxy_catalog(
        cluster_m[i],
        cluster_z[i],
        concentration[i],
        cosmo,
        cluster_ra=ra[i],
        cluster_dec=dec[i],
        zsrc="chang13",
        delta_so=200,
        massdef="critical",
        halo_profile_model="nfw",
        zsrc_min=cluster_z[i] + 0.1,
        zsrc_max=3.0,
        field_size=10.0,
        shapenoise=0.04,
        photoz_sigma_unscaled=0.02,
        ngals=n_gals,
        mean_e_err=0.1,
    )

    cl = clmm.GalaxyCluster("mock_cluster", ra[i], dec[i], cluster_z[i], noisy_data_z)

    # compute DeltaSigma for each background galaxy
    cl.compute_tangential_and_cross_components(
        shape_component1="e1",
        shape_component2="e2",
        tan_component="DS_t",
        cross_component="DS_x",
        cosmo=cosmo,
        is_deltasigma=True,
        use_pdz=True,
    )

    # compute the weights to be used to bluid the DeltaSigma radial profiles
    cl.compute_galaxy_weights(
        use_pdz=True,
        use_shape_noise=True,
        shape_component1="e1",
        shape_component2="e2",
        use_shape_error=True,
        shape_component1_err="e_err",
        shape_component2_err="e_err",
        weight_name="w_ls",
        cosmo=cosmo,
        is_deltasigma=True,
        add=True,
    )

    # append the cluster in the list
    gclist.append(cl)

cl.galcat.columns

<TableColumns names=('ra','dec','e1','e2','e_err','z','ztrue','pzpdf','id','sigma_c_eff','theta','DS_t','DS_x','w_ls')>

Creating ClusterEnsemble object and estimation of individual excess surface density profiles

From the galaxy cluster object list gclist, we instantiate a cluster ensemble object clusterensemble. This instantiation step uses - the individual galaxy input \(\Delta\Sigma_+\) and \(\Delta\Sigma_{\times}\) values (computed in the previous step, DS_{t,x}) - the corresponding individual input weights w_ls (computed in the previous step)

to compute

  • the output tangential DS_t and cross signal DS_x binned profiles (where the binning is controlled by bins)

  • the associated profile weights W_l (that will be used to compute the stacked profile)

ensemble_id = 1
names = ["id", "ra", "dec", "z", "radius", "gt", "gx", "W_l"]
bins = np.logspace(np.log10(0.3), np.log10(5), 10)
clusterensemble = ClusterEnsemble(
    ensemble_id,
    gclist,
    tan_component_in="DS_t",
    cross_component_in="DS_x",
    tan_component_out="DS_t",
    cross_component_out="DS_x",
    weights_in="w_ls",
    weights_out="W_l",
    bins=bins,
    bin_units="Mpc",
    cosmo=cosmo,
)

There is also the option to create an ClusterEnsemble object without the clusters list. Then, the user may compute the individual profile for each wanted cluster and compute the radial profile once all the indvidual profiles have been computed. This method may be reccomended if there a large number of clusters to avoid excess of memory allocation, since the user can generate each cluster separately, compute its individual profile and then delete the cluster object.

ensemble_id2 = 2
empty_cluster_ensemble = ClusterEnsemble(ensemble_id2)
for cluster in gclist:
    empty_cluster_ensemble.make_individual_radial_profile(
        galaxycluster=cluster,
        tan_component_in="DS_t",
        cross_component_in="DS_x",
        tan_component_out="DS_t",
        cross_component_out="DS_x",
        weights_in="w_ls",
        weights_out="W_l",
        bins=bins,
        bin_units="Mpc",
        cosmo=cosmo,
    )

A third option is where all clusters already have the profile computed, and we just add those to the ClusterEnsemble object:

# add profiles to gclist
for cluster in gclist:
    cluster.make_radial_profile(
        tan_component_in="DS_t",
        cross_component_in="DS_x",
        tan_component_out="DS_t",
        cross_component_out="DS_x",
        weights_in="w_ls",
        weights_out="W_l",
        bins=bins,
        bin_units="Mpc",
        cosmo=cosmo,
        use_weights=True,
    )

cluster.profile
GCData
defined by: cosmo='CCLCosmology(H0=71.0, Omega_dm0=0.2202, Omega_b0=0.0448, Omega_k0=0.0)', bin_units='Mpc'
with columns: radius_min, radius, radius_max, DS_t, DS_t_err, DS_x, DS_x_err, z, z_err, n_src, W_l
9 objects
radius_minradiusradius_maxDS_tDS_t_errDS_xDS_x_errzz_errn_srcW_l
float64float64float64float64float64float64float64float64float64int64float64
0.299999999999999930.35538888084859040.410092877858960247114109394812.8974981653954774.0371725831060887.9770821817988064.01.22941885527082380.1254222511654474296.36216161623601e-29
0.41009287785896020.49712754212398570.560587228235480428407616591812.97742764109085226.86-57071960817478.7246736076941033.971.3771445318850360.09702952618749792491.15161521880834e-28
0.56058722823548040.67123612279400550.766309432393553274964576073867.1938525941380754.6456308907913134.0242228781968987.931.40053883231784870.06709616391264094862.113338478541802e-28
0.76630943239355320.90638180342879091.047526801535795231222548973785.8628953248329270.2034403827566750.54330962934549207.2421.3571852747152190.0547948614357719951583.639913884168516e-28
1.04752680153579521.25737626108032761.431944268920687443299578835108.5421990385325328.957150722770699.8338620759010434230.5271.34516410410421150.034977356940368692967.038287482539166e-28
1.43194426892068741.70438577640466861.95743382058443269749022992700.26415940514809587.195-8821468381838.70115585239988553.3791.33063787396155030.026819503616418075351.270287908236562e-27
1.95743382058443262.3390498533169922.675765562339765626294486206933.4710815129177769.555-29746134079554.89510982086708852.971.35951921466047550.0195703169606892610512.521188293741188e-27
2.67576556233976563.1943994280425953.65770799978604534645975281894.5787782221363303.094-6755554338437.7027744746597065.1461.36790353652810180.01453405310911619319514.719527212638533e-27
3.65770799978604534.3684553488155725.0000000000000013195312396468.36675740366240659.977-1692230857239.3755890154465299.2721.34960897269319570.01051691211193225536448.730103338179242e-27
ensemble_id3 = 3
empty_cluster_ensemble = ClusterEnsemble(ensemble_id3)
for cluster in gclist:
    empty_cluster_ensemble.add_individual_radial_profile(
        galaxycluster=cluster,
        profile_table=cluster.profile,
        tan_component="DS_t",
        cross_component="DS_x",
        weights="W_l",
    )