Boost factors
=============

.. code:: ipython3

    import clmm
    from clmm import Cosmology
    from clmm.support import mock_data as mock
    from clmm.galaxycluster import GalaxyCluster
    import matplotlib.pyplot as plt
    import sys
    import clmm.utils as u

Make sure we know which version we’re using

.. code:: ipython3

    clmm.__version__




.. parsed-literal::

    '1.12.0'



Define cosmology object
~~~~~~~~~~~~~~~~~~~~~~~

.. code:: ipython3

    mock_cosmo = Cosmology(H0=70.0, Omega_dm0=0.27 - 0.045, Omega_b0=0.045, Omega_k0=0.0)

First, we want to generate a :math:`\Delta\Sigma` (excess surface
density) profile from mock data, to which we can apply boost factors.
The mock data is generated in the following cells.

Generate cluster object from mock data

.. code:: ipython3

    cosmo = mock_cosmo
    cluster_id = "Awesome_cluster"
    cluster_m = 1.0e15
    cluster_z = 0.3
    concentration = 4
    ngals = 1000
    
    zsrc_min = cluster_z + 0.1  # we only want to draw background galaxies
    
    noisy_data_z = mock.generate_galaxy_catalog(
        cluster_m,
        cluster_z,
        concentration,
        cosmo,
        "chang13",
        zsrc_min=zsrc_min,
        shapenoise=0.005,
        photoz_sigma_unscaled=0.05,
        ngals=ngals,
    )

Loading this into a CLMM cluster object centered on (0,0)

.. code:: ipython3

    cluster_ra = 0.0
    cluster_dec = 0.0
    cl = GalaxyCluster(cluster_id, cluster_ra, cluster_dec, cluster_z, noisy_data_z)

Compute cross and tangential excess surface density for each source
galaxy

.. code:: ipython3

    _ = cl.compute_tangential_and_cross_components(
        geometry="flat",
        shape_component1="e1",
        shape_component2="e2",
        tan_component="DeltaSigma_tan",
        cross_component="DeltaSigma_cross",
        add=True,
        cosmo=cosmo,
        is_deltasigma=True,
    )


.. parsed-literal::

    /global/homes/a/aguena/.local/perlmutter/python-3.11/lib/python3.11/site-packages/clmm-1.12.0-py3.11.egg/clmm/cosmology/parent_class.py:110: UserWarning: 
    Some source redshifts are lower than the cluster redshift.
    Sigma_crit = np.inf for those galaxies.


Calculate the binned profile

.. code:: ipython3

    cl.make_radial_profile(
        "kpc",
        cosmo=cosmo,
        tan_component_in="DeltaSigma_tan",
        cross_component_in="DeltaSigma_cross",
        tan_component_out="DeltaSigma_tan",
        cross_component_out="DeltaSigma_cross",
        table_name="DeltaSigma_profile",
    )
    # Format columns for display
    for col in cl.DeltaSigma_profile.colnames:
        fmt = cl.DeltaSigma_profile[col].info.format
        if "DeltaSigma" in col:
            fmt = ".2e"
        elif any(typ in col for typ in ("z", "radius")):
            fmt = ".2f"
        cl.DeltaSigma_profile[col].info.format = fmt
    # Show
    cl.DeltaSigma_profile.show_in_notebook()




.. raw:: html

    <i>GCData length=10</i>
    <table id="table139955615623888-471625" class="table-striped table-bordered table-condensed">
    <thead><tr><th>idx</th><th>radius_min</th><th>radius</th><th>radius_max</th><th>DeltaSigma_tan</th><th>DeltaSigma_tan_err</th><th>DeltaSigma_cross</th><th>DeltaSigma_cross_err</th><th>z</th><th>z_err</th><th>n_src</th><th>W_l</th></tr></thead>
    <tr><td>0</td><td>116.87</td><td>441.96</td><td>652.82</td><td>2.38e+14</td><td>3.64e+13</td><td>-1.24e+13</td><td>1.62e+13</td><td>1.19</td><td>0.14</td><td>18</td><td>18.0</td></tr>
    <tr><td>1</td><td>652.82</td><td>967.93</td><td>1188.77</td><td>1.16e+14</td><td>3.55e+12</td><td>1.56e+11</td><td>2.06e+12</td><td>1.17</td><td>0.10</td><td>47</td><td>47.0</td></tr>
    <tr><td>2</td><td>1188.77</td><td>1480.84</td><td>1724.72</td><td>7.26e+13</td><td>2.05e+12</td><td>2.97e+12</td><td>2.51e+12</td><td>1.37</td><td>0.09</td><td>67</td><td>67.0</td></tr>
    <tr><td>3</td><td>1724.72</td><td>1989.06</td><td>2260.67</td><td>5.70e+13</td><td>2.01e+12</td><td>-3.14e+11</td><td>1.72e+12</td><td>1.30</td><td>0.08</td><td>92</td><td>92.0</td></tr>
    <tr><td>4</td><td>2260.67</td><td>2531.03</td><td>2796.62</td><td>4.57e+13</td><td>1.80e+12</td><td>-1.45e+12</td><td>1.38e+12</td><td>1.41</td><td>0.07</td><td>139</td><td>139.0</td></tr>
    <tr><td>5</td><td>2796.62</td><td>3077.13</td><td>3332.57</td><td>3.94e+13</td><td>2.97e+12</td><td>3.13e+11</td><td>1.77e+12</td><td>1.27</td><td>0.06</td><td>167</td><td>167.0</td></tr>
    <tr><td>6</td><td>3332.57</td><td>3598.82</td><td>3868.52</td><td>3.03e+13</td><td>2.02e+12</td><td>2.58e+12</td><td>3.69e+12</td><td>1.21</td><td>0.05</td><td>196</td><td>196.0</td></tr>
    <tr><td>7</td><td>3868.52</td><td>4127.03</td><td>4404.48</td><td>2.58e+13</td><td>1.56e+12</td><td>1.54e+12</td><td>1.26e+12</td><td>1.21</td><td>0.05</td><td>163</td><td>163.0</td></tr>
    <tr><td>8</td><td>4404.48</td><td>4659.49</td><td>4940.43</td><td>2.26e+13</td><td>1.80e+12</td><td>3.08e+12</td><td>1.74e+12</td><td>1.37</td><td>0.09</td><td>76</td><td>76.0</td></tr>
    <tr><td>9</td><td>4940.43</td><td>5184.62</td><td>5476.38</td><td>2.37e+13</td><td>9.03e+12</td><td>-2.74e+11</td><td>3.98e+12</td><td>1.14</td><td>0.13</td><td>35</td><td>35.0</td></tr>
    </table><style>table.dataTable {clear: both; width: auto !important; margin: 0 !important;}
    .dataTables_info, .dataTables_length, .dataTables_filter, .dataTables_paginate{
    display: inline-block; margin-right: 1em; }
    .paginate_button { margin-right: 5px; }
    </style>
    <script>
    
    var astropy_sort_num = function(a, b) {
        var a_num = parseFloat(a);
        var b_num = parseFloat(b);
    
        if (isNaN(a_num) && isNaN(b_num))
            return ((a < b) ? -1 : ((a > b) ? 1 : 0));
        else if (!isNaN(a_num) && !isNaN(b_num))
            return ((a_num < b_num) ? -1 : ((a_num > b_num) ? 1 : 0));
        else
            return isNaN(a_num) ? -1 : 1;
    }
    
    require.config({paths: {
        datatables: 'https://cdn.datatables.net/1.10.12/js/jquery.dataTables.min'
    }});
    require(["datatables"], function(){
        console.log("$('#table139955615623888-471625').dataTable()");
    
    jQuery.extend( jQuery.fn.dataTableExt.oSort, {
        "optionalnum-asc": astropy_sort_num,
        "optionalnum-desc": function (a,b) { return -astropy_sort_num(a, b); }
    });
    
        $('#table139955615623888-471625').dataTable({
            order: [],
            pageLength: 50,
            lengthMenu: [[10, 25, 50, 100, 500, 1000, -1], [10, 25, 50, 100, 500, 1000, 'All']],
            pagingType: "full_numbers",
            columnDefs: [{targets: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11], type: "optionalnum"}]
        });
    });
    </script>




Plot the :math:`\Delta\Sigma` profile

.. code:: ipython3

    plt.errorbar(
        cl.DeltaSigma_profile["radius"],
        cl.DeltaSigma_profile["DeltaSigma_tan"],
        cl.DeltaSigma_profile["DeltaSigma_tan_err"],
        marker="o",
    )
    
    plt.title("DeltaSigma profile")
    plt.xlabel("Radius [kpc]")
    plt.ylabel("$\Delta\Sigma [M_\odot\; Mpc^{-2}]$")
    
    plt.show()



.. image:: demo_boost_factors_files/demo_boost_factors_16_0.png


Boost Factors
-------------

CLMM offers two boost models, the NFW boost model, and a powerlaw boost
model.

Note that ``compute_nfw_boost`` requires two parameters to be specified,
``rs`` and ``b0``, and ``compute_powerlaw_boost`` requires three
paramters, ``rs``, ``b0`` and ``alpha``. The default values are in kpc.

Details on these boost models can be found
`here <https://cluster-toolkit.readthedocs.io/en/latest/source/boostfactors.html>`__

First, we can calculate the boost factors from the two models.

.. code:: ipython3

    nfw_boost = u.compute_nfw_boost(cl.DeltaSigma_profile["radius"], rscale=1000, boost0=0.1)
    
    powerlaw_boost = u.compute_powerlaw_boost(
        cl.DeltaSigma_profile["radius"], rscale=1000, boost0=0.1, alpha=-1.0
    )

Plot the two boost factors, :math:`B(R)`

.. code:: ipython3

    plt.plot(cl.DeltaSigma_profile["radius"], nfw_boost, label="NFW boost factor")
    plt.plot(cl.DeltaSigma_profile["radius"], powerlaw_boost, label="Powerlaw boost factor")
    plt.xlabel("Radius [kpc]")
    plt.ylabel("$B(R)$")
    plt.legend()
    plt.show()



.. image:: demo_boost_factors_files/demo_boost_factors_21_0.png


The :math:`\Delta\Sigma` profiles can be corrected with the boost factor
using ``correct_sigma_with_boost_values`` or
``correct_sigma_with_boost_model``.

``correct_sigma_with_boost_values`` requires us to precompute the boost
factor, e.g. using ``compute_nfw_boost``.
``correct_sigma_with_boost_model`` simply requires us to specify the
boost model.

Note that the boost factor can be used in one of two ways.

Either the boost factor can be applied to the observed data vector to
correct for the dilution of the signal by cluster member galaxies. In
this case the amplitude of the corrected profile will increase.

Or the boost factor can be applied to the model prediction. In this case
it behaves as a dilution factor, and the resulting model prediction will
be lower than the original one.

Both scenarios will improve the agreement between the mock data and
observed data, by accounting for cluster member galaxy contamination.

In this notebook, we use the second approach, where the data is
generated using mock data that does not account for dilution until the
boost factor is applied. The corrected profiles from the mock data are
lower than the uncorrected one.

Essentially we are diluting the mock profile to mimick the effect of
contamination by cluster members.

First we will apply the boost factor with
``correct_sigma_with_boost_values``

.. code:: ipython3

    Sigma_corrected_powerlaw_boost = u.correct_sigma_with_boost_values(
        cl.DeltaSigma_profile["DeltaSigma_tan"], powerlaw_boost
    )
    Sigma_corrected_nfw_boost = u.correct_sigma_with_boost_values(
        cl.DeltaSigma_profile["DeltaSigma_tan"], nfw_boost
    )

Plot the result

.. code:: ipython3

    plt.errorbar(
        cl.DeltaSigma_profile["radius"],
        Sigma_corrected_nfw_boost,
        cl.DeltaSigma_profile["DeltaSigma_tan_err"],
        marker="o",
        label="$\Delta \Sigma$ / NFW boost factor",
    )
    
    plt.errorbar(
        cl.DeltaSigma_profile["radius"],
        Sigma_corrected_powerlaw_boost,
        cl.DeltaSigma_profile["DeltaSigma_tan_err"],
        marker="o",
        label="$\Delta \Sigma$ / Powerlaw boost factor",
    )
    
    plt.errorbar(
        cl.DeltaSigma_profile["radius"],
        cl.DeltaSigma_profile["DeltaSigma_tan"],
        cl.DeltaSigma_profile["DeltaSigma_tan_err"],
        marker="o",
        label="uncorrected $\Delta \Sigma$",
        color="grey",
    )
    
    # plt.loglog()
    plt.title("DeltaSigma profile")
    plt.xlabel("Radius [kpc]")
    plt.ylabel("$\Delta\Sigma [M_\odot\; Mpc^{-2}]$")
    plt.legend()
    plt.show()



.. image:: demo_boost_factors_files/demo_boost_factors_27_0.png


Now the same again but with ``correct_sigma_with_boost_model``

.. code:: ipython3

    Sigma_corrected_powerlaw_boost = u.correct_sigma_with_boost_model(
        cl.DeltaSigma_profile["radius"],
        cl.DeltaSigma_profile["DeltaSigma_tan"],
        boost_model="powerlaw_boost",
    )
    Sigma_corrected_nfw_boost = u.correct_sigma_with_boost_model(
        cl.DeltaSigma_profile["radius"],
        cl.DeltaSigma_profile["DeltaSigma_tan"],
        boost_model="nfw_boost",
    )
    
    plt.errorbar(
        cl.DeltaSigma_profile["radius"],
        Sigma_corrected_nfw_boost,
        cl.DeltaSigma_profile["DeltaSigma_tan_err"],
        marker="o",
        label="$\Delta \Sigma$ / NFW boost factor",
    )
    
    plt.errorbar(
        cl.DeltaSigma_profile["radius"],
        Sigma_corrected_powerlaw_boost,
        cl.DeltaSigma_profile["DeltaSigma_tan_err"],
        marker="o",
        label="$\Delta \Sigma$ / Powerlaw boost factor",
    )
    
    plt.errorbar(
        cl.DeltaSigma_profile["radius"],
        cl.DeltaSigma_profile["DeltaSigma_tan"],
        cl.DeltaSigma_profile["DeltaSigma_tan_err"],
        marker="o",
        label="uncorrected $\Delta \Sigma$",
        color="grey",
    )
    
    # plt.loglog()
    plt.title("DeltaSigma profile")
    plt.xlabel("Radius [kpc]")
    plt.ylabel("$\Delta\Sigma [M_\odot\; Mpc^{-2}]$")
    plt.legend()
    plt.show()



.. image:: demo_boost_factors_files/demo_boost_factors_29_0.png