theory_utils

multiprobe_framework.theory_utils

get_lss_spectra(cosmo_ccl, theta_dict, spectra, cl_dict, ells, path_boss, path_kids, path_desi=None)

Setup the parameters for the CCL calculations.

Note that the clustering number counts tracers currently do not include any magnification bias as these are assumed to be only large for ell<20. This is a simplification that can be changed in the future.

Parameters:

Name	Type	Description	Default
theta_dict		Dictionary containing the cosmological parameters	required
path_boss		Path to the BOSS data	required
path_kids		Path to the KIDS data	required
spectra		List of spectra to calculate	required

Returns:

Type	Description
	cl_params_list: List of dictionaries containing the parameters for the pycosmo calculations

Source code in src/multiprobe_framework/theory_utils.py

def get_lss_spectra(
    cosmo_ccl,
    theta_dict,
    spectra,
    cl_dict,
    ells,
    path_boss,
    path_kids,
    path_desi=None,  # Path to DESI data if adding this in
):
    """
    Setup the parameters for the CCL calculations.

    Note that the clustering number counts tracers currently do not include
    any magnification bias as these are assumed to be only large for ell<20.
    This is a simplification that can be changed in the future.

    :param theta_dict: Dictionary containing the cosmological parameters
    :param path_boss: Path to the BOSS data
    :param path_kids: Path to the KIDS data
    :param spectra: List of spectra to calculate

    :return: cl_params_list: List of dictionaries containing the parameters for
                             the pycosmo calculations
    """
    # Paths to the redshift distributions
    nz_lowz = np.loadtxt(f"{path_boss}/nz_lowz_0p0_0p4.txt")
    nz_cmass = np.loadtxt(f"{path_boss}/nz_cmass_0p4_0p8.txt")
    kids_redshift_dict = {
        i: np.loadtxt(f"{path_kids}/nz_{i}_no_err") for i in range(1, 6)
    }

    # delete any below 0 after del_z shifts and make all lensing redshift
    # distributions start at 1e-6 for numerical stability
    z_kids = {}
    nz_kids = {}
    for i in range(1, 6):
        z = kids_redshift_dict[i][:, 0] + theta_dict[f"del_z{i}"]
        mask = z > 0
        z_kids[i] = z[mask]
        nz_kids[i] = kids_redshift_dict[i][:, 1][mask]
        if z_kids[i][0] < 1e-6:
            z_kids[i][0] = 1e-6

    if path_desi is not None:
        nz_desi_dict = {i: np.loadtxt(f"{path_desi}/nz_{i}.txt") for i in range(1, 5)}
        z_desi = {}
        nz_desi = {}
        for i in range(1, 5):
            z = nz_desi_dict[i][:, 0]
            nz = nz_desi_dict[i][:, 1]
            z_desi[i] = z
            nz_desi[i] = nz

        # from Table 1. of NS paper
        mag_default_dict = {"d1": 0.972, "d2": 1.044, "d3": 0.974, "d4": 0.988}

        if "mag_bias_d1" in theta_dict:
            mag_dict = {f"d{i}": theta_dict[f"mag_bias_d{i}"] for i in range(1, 5)}

        else:
            mag_dict = {f"d{i}": mag_default_dict[f"d{i}"] for i in range(1, 5)}

    tracer_dict = {
        "l": ccl.NumberCountsTracer(
            cosmo_ccl,
            has_rsd=True,
            dndz=(nz_lowz[:, 0], nz_lowz[:, 1]),
            bias=(nz_lowz[:, 0], theta_dict["bl"] * np.ones_like(nz_lowz[:, 0])),
        ),
        "c": ccl.NumberCountsTracer(
            cosmo_ccl,
            has_rsd=True,
            dndz=(nz_cmass[:, 0], nz_cmass[:, 1]),
            bias=(nz_lowz[:, 0], theta_dict["bc"] * np.ones_like(nz_lowz[:, 0])),
        ),
        "d1": ccl.NumberCountsTracer(
            cosmo_ccl,
            has_rsd=True,
            dndz=(z_desi[1], nz_desi[1]),
            bias=(z_desi[1], theta_dict["bias_d1"] * np.ones_like(z_desi[1])),
            mag_bias=(z_desi[1], mag_dict["d1"] * np.ones_like(z_desi[1])),
        ),
        "d2": ccl.NumberCountsTracer(
            cosmo_ccl,
            has_rsd=True,
            dndz=(z_desi[2], nz_desi[2]),
            bias=(z_desi[2], theta_dict["bias_d2"] * np.ones_like(z_desi[2])),
            mag_bias=(z_desi[2], mag_dict["d2"] * np.ones_like(z_desi[1])),
        ),
        "d3": ccl.NumberCountsTracer(
            cosmo_ccl,
            has_rsd=True,
            dndz=(z_desi[3], nz_desi[3]),
            bias=(z_desi[3], theta_dict["bias_d3"] * np.ones_like(z_desi[3])),
            mag_bias=(z_desi[3], mag_dict["d3"] * np.ones_like(z_desi[1])),
        ),
        "d4": ccl.NumberCountsTracer(
            cosmo_ccl,
            has_rsd=True,
            dndz=(z_desi[4], nz_desi[4]),
            bias=(z_desi[4], theta_dict["bias_d4"] * np.ones_like(z_desi[4])),
            mag_bias=(z_desi[4], mag_dict["d4"] * np.ones_like(z_desi[1])),
        ),
        "g1": ccl.WeakLensingTracer(
            cosmo_ccl,
            dndz=(z_kids[1], nz_kids[1]),
            use_A_ia=True,
            ia_bias=(z_kids[1], theta_dict["A_IA"] * np.ones_like(nz_kids[1])),
        ),
        "g2": ccl.WeakLensingTracer(
            cosmo_ccl,
            dndz=(z_kids[2], nz_kids[2]),
            use_A_ia=True,
            ia_bias=(z_kids[2], theta_dict["A_IA"] * np.ones_like(nz_kids[2])),
        ),
        "g3": ccl.WeakLensingTracer(
            cosmo_ccl,
            dndz=(z_kids[3], nz_kids[3]),
            use_A_ia=True,
            ia_bias=(z_kids[3], theta_dict["A_IA"] * np.ones_like(nz_kids[3])),
        ),
        "g4": ccl.WeakLensingTracer(
            cosmo_ccl,
            dndz=(z_kids[4], nz_kids[4]),
            use_A_ia=True,
            ia_bias=(z_kids[4], theta_dict["A_IA"] * np.ones_like(nz_kids[4])),
        ),
        "g5": ccl.WeakLensingTracer(
            cosmo_ccl,
            dndz=(z_kids[5], nz_kids[5]),
            use_A_ia=True,
            ia_bias=(z_kids[5], theta_dict["A_IA"] * np.ones_like(nz_kids[5])),
        ),
        "k": ccl.CMBLensingTracer(
            cosmo_ccl, z_source=1090.0, n_samples=1000
        ),  # increase n_samples to 5X default for good accuracy
        "ka": ccl.CMBLensingTracer(
            cosmo_ccl, z_source=1090.0, n_samples=1000
        ),  # increase n_samples to 5X default for good accuracy
        "t": ccl.ISWTracer(cosmo_ccl, n_chi=4096),
    }

    for spec in spectra:
        if spec not in ["tt", "te", "ee"]:
            specind_1, specind_2 = utils.get_specinds(spec)
            tracer_1 = tracer_dict[specind_1]
            tracer_2 = tracer_dict[specind_2]

            # Check if either tracer is ISW (tracer key 't')
            # For ISW cross-correlations, use standard Limber (higher l_limber value)
            if specind_1 == "t" or specind_2 == "t":
                l_limber_value = (
                    -1
                )  # Use standard Limber for all ell (no beyond-Limber)
            else:
                l_limber_value = (
                    60  # Use beyond Limber for ell < 70, Limber for ell >= 70
                )

            cl_dict[spec] = ccl.angular_cl(
                cosmo_ccl,
                tracer_1,
                tracer_2,
                ells,
                l_limber=l_limber_value,
                limber_integration_method="qag_quad",  # More accurate than 'spline'
                non_limber_integration_method="FKEM",  # Control beyond-Limber accuracy
            )

    return cl_dict

apply_binning(spec, path_to_workspace, cl_dict)

Reads and processes the spectrum based on the spec identifier.

Parameters:

Name	Type	Description	Default
spec		The spectrum identifier	required
path_to_workspace		Path to the NaMaster workspace	required
cl_dict		Dictionary containing the fullsky cls	required

Returns:

Type	Description
	cl_realisation: The binned cls

Source code in src/multiprobe_framework/theory_utils.py

def apply_binning(spec, path_to_workspace, cl_dict):
    """
    Reads and processes the spectrum based on the spec identifier.

    :param spec: The spectrum identifier
    :param path_to_workspace: Path to the NaMaster workspace
    :param cl_dict: Dictionary containing the fullsky cls

    :return: cl_realisation: The binned cls
    """
    if spec in ["tt", "ee", "te", "kk", "kaka"]:
        return cl_dict[f"{spec}"]  # Return cl_full only, no binning of cmb cls

    w = nmt.NmtWorkspace()
    w.read_from(f"{path_to_workspace}/pymaster_workspace_{spec}.fits")

    cl_full = cl_dict[spec]

    # Set monopole and dipole to zero for consistency with rest of analysis.
    cl_full[:2] = 0

    # Configure the cell coupling depending on the spectrum indices
    specind_1, specind_2 = utils.get_specinds(spec)

    def is_shear(index):
        return index in ["g1", "g2", "g3", "g4", "g5"]

    if is_shear(specind_1) and is_shear(specind_2):
        coupled = w.couple_cell(
            np.array(
                [
                    np.array(cl_full),
                    np.zeros_like(cl_full),
                    np.zeros_like(cl_full),
                    np.zeros_like(cl_full),
                ]
            )
        )
    elif not is_shear(specind_1) and not is_shear(specind_2):
        coupled = w.couple_cell(np.array([cl_full]))
    else:
        coupled = w.couple_cell([np.array(cl_full), np.zeros_like(cl_full)])

    cl_realisation = w.decouple_cell(coupled)
    return cl_realisation[0]