diff --git a/sprint_week_2026/DP2_cluster_mass_fit_color_selection-cleaned.ipynb b/sprint_week_2026/DP2_cluster_mass_fit_color_selection-cleaned.ipynb new file mode 100644 index 0000000..beb13f7 --- /dev/null +++ b/sprint_week_2026/DP2_cluster_mass_fit_color_selection-cleaned.ipynb @@ -0,0 +1,890 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Galaxy cluster weak-lensing mass with DP2 and CLMM\n", + "\n", + "For the Rubin Science Platform at data.lsst.cloud.
\n", + "Data Release: Data Preview 2
\n", + "Container Size: large
\n", + "LSST Science Pipelines version: r30 (same as DP2 tutorial 304.3)
\n", + "Packages: `lsst.rsp`, `clmm` (≥ 1.18), a CLMM theory backend (`pyccl`)\n", + "\n", + "**Learning objective:** Measure the reduced tangential shear profile of the galaxy cluster\n", + "PSZ2 G309.43-72.86 with DP2 data and fit its halo mass with CLMM.\n", + "\n", + "**Credit:** Data access and galaxy selection follow the Rubin DP2 tutorial\n", + "[304.3 \"Galaxy cluster weak lensing with HSM\"](https://dp2.lsst.io/tutorials/notebook/304/notebook-304-3.html).\n", + "The mass-fitting approach follows the CLMM example\n", + "`examples/mass_fitting/Example3_Fit_Halo_Mass_to_Shear_Catalog.ipynb`." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Overview\n", + "\n", + "PSZ2 G309.43-72.86 is a rich galaxy cluster at redshift $z=0.35$ in DP2.\n", + "Tutorial 304.3 shows that background galaxies around it have a positive mean *tangential\n", + "ellipticity*. Here we go one step further and turn that signal into a **cluster mass**:\n", + "\n", + "1. **Setup** – make CLMM importable on the RSP.\n", + "2. **Load the DP2 Object catalog** around the cluster with TAP (as in tutorial 304.3).\n", + " This is the only slow step, and it is run only once.\n", + "3. **Define the analysis parameters** in one place.\n", + "4. **Select galaxies** (extended, well-measured objects).\n", + "5. **Select background galaxies by color** (redder than the cluster red sequence).\n", + "6. **Apply shape-quality cuts.**\n", + "7. **Convert HSM ellipticities into shear estimates.**\n", + "8. **Compute tangential and cross shear** with CLMM.\n", + "9. **Build the radial shear profile** in physical units (Mpc).\n", + "10. **Assume a redshift distribution** for the background sources (no photo-z).\n", + "11. **Fit an NFW halo mass** with CLMM.\n", + "12. **Plot the result.**\n", + "13. **Caveats and exercises.**\n", + "\n", + "**Tip:** the catalog from Section 2 is never modified afterwards. To try other choices (for example\n", + "the radial range), change the parameters in Section 3 and rerun from Section 3 onwards, without\n", + "querying again." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1. Setup\n", + "\n", + "### 1.1 One-time installation (only the first time)\n", + "\n", + "CLMM is not part of the LSST Science Pipelines, so we use the latest CLMM code from GitHub.\n", + "Open a terminal on the RSP, go to any directory you like, and clone CLMM:\n", + "\n", + "```bash\n", + "git clone https://github.com/LSSTDESC/CLMM.git\n", + "```\n", + "\n", + "To get newer CLMM code later, run `git pull` inside the `CLMM` folder.\n", + "\n", + "`import clmm` also needs `qp` (used by CLMM for photo-z PDFs; imported even though this notebook\n", + "does not use photo-z) and a theory backend (we use CCL). If they are not available in your\n", + "RSP environment, uncomment and run the cell below **once**, then restart the kernel\n", + "(Kernel → Restart Kernel)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# %pip install --user qp-prob pyccl" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 1.2 Make CLMM importable\n", + "\n", + "Jupyter runs a notebook from the notebook's own folder, and the CLMM clone's folder needs to be provided. We put it at the front of the Python path, so `import clmm` uses this\n", + "GitHub version instead of any older installed version. We also ask CLMM to use the CCL backend." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import sys\n", + "\n", + "CLMM_PATH = os.path.abspath(\"./CLMM\") # The clone's top folder, which contains the \"clmm\" package\n", + "sys.path.insert(0, CLMM_PATH)\n", + "os.environ[\"CLMM_MODELING_BACKEND\"] = \"ccl\" # Must be set before importing clmm\n", + "\n", + "import clmm\n", + "\n", + "print(\"CLMM version:\", clmm.__version__)\n", + "print(\"CLMM loaded from:\", clmm.__file__)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 1.3 Import packages\n", + "\n", + "General packages (numpy, scipy, matplotlib, astropy), the RSP data-access package (`lsst.rsp`)\n", + "and the CLMM modules we need." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from astropy.coordinates import SkyCoord\n", + "from lsst.rsp import RSPDiscovery\n", + "from scipy.integrate import quad\n", + "\n", + "import clmm.dataops as da\n", + "from clmm.support.sampler import fitters\n", + "from clmm.utils import (\n", + " compute_beta_s_mean_from_distribution,\n", + " compute_beta_s_square_mean_from_distribution,\n", + " make_bins,\n", + " redshift_distributions,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Make the axis labels and tick labels in all figures 1.7 times larger than the matplotlib default." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "label_fontsize = 1.7 * plt.rcParams[\"font.size\"]\n", + "plt.rcParams.update(\n", + " {\n", + " \"axes.labelsize\": label_fontsize,\n", + " \"xtick.labelsize\": label_fontsize,\n", + " \"ytick.labelsize\": label_fontsize,\n", + " }\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2. Load the DP2 Object catalog\n", + "\n", + "This section is run **once**. The query only depends on the cluster position (we use the\n", + "Brightest Cluster Galaxy, BCG, as the center), the band used for galaxy shapes, and the search\n", + "radius." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Cluster PSZ2 G309.43-72.86 (BCG position used as the center)\n", + "cluster_ra = 10.2082 # [deg]\n", + "cluster_dec = -44.1307 # [deg]\n", + "\n", + "band_shape = \"i\" # Band used for galaxy shapes\n", + "radius_deg = 0.4 # Cone-search radius around the BCG [deg]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Get the TAP service, exactly as in tutorial 304.3." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "discovery = RSPDiscovery(\"dp2\")\n", + "service = discovery.get_tap_client()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Query the `dp2.Object` table for objects within `radius_deg` of the BCG. We ask for the HSM\n", + "Regaussianization shapes in the `band_shape` band, the riz CModel photometry, and quality flags.\n", + "The signal-to-noise ratio (SNR) is computed in the query because it is not a column in the schema." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "query = (\n", + " \"SELECT objectId, coord_ra, coord_dec, refExtendedness, \"\n", + " f\"{band_shape}_blendedness, \"\n", + " f\"{band_shape}_hsmShapeRegauss_e1, \"\n", + " f\"{band_shape}_hsmShapeRegauss_e2, \"\n", + " f\"{band_shape}_hsmShapeRegauss_sigma, \"\n", + " f\"{band_shape}_hsmShapeRegauss_flag, \"\n", + " \"r_cModelMag, i_cModelMag, z_cModelMag, \"\n", + " \"r_cModelFlux/r_cModelFluxErr AS r_cModel_SNR, \"\n", + " \"i_cModelFlux/i_cModelFluxErr AS i_cModel_SNR, \"\n", + " \"z_cModelFlux/z_cModelFluxErr AS z_cModel_SNR, \"\n", + " \"r_cModel_flag, i_cModel_flag, z_cModel_flag \"\n", + " \"FROM dp2.Object \"\n", + " \"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), \"\n", + " f\"CIRCLE('ICRS', {cluster_ra}, {cluster_dec}, {radius_deg})) = 1 \"\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Run the query, check that it completed, and convert the result into an astropy table." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "job = service.submit_job(query)\n", + "job.run()\n", + "job.wait(phases=[\"COMPLETED\", \"ERROR\"])\n", + "print(\"Job phase is\", job.phase)\n", + "if job.phase == \"ERROR\":\n", + " job.raise_if_error()\n", + "assert job.phase == \"COMPLETED\"\n", + "obj_cat = job.fetch_result().to_table()\n", + "print(\"Number of objects is\", len(obj_cat))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 3. Define the analysis parameters\n", + "\n", + "All analysis choices are collected here. **After changing any of them, rerun from this cell\n", + "onwards**; there is no need to query the catalog again.\n", + "\n", + "* **Cluster redshift:** $z = 0.35$.\n", + "* **Selection:** the same values as tutorial 304.3.\n", + "* **Shear calibration:** DP2 does not provide multiplicative ($m$) or additive ($c$) shear biases\n", + " for the HSM shapes, so we set $m = c = 0$. We assume a per-component shape dispersion\n", + " $e_{\\rm rms} = 0.365$.\n", + "* **Radial bins:** the same bins as tutorial 304.3, which uses 5 logarithmic bins from 4 to 24 arcmin.\n", + " At $z = 0.35$ in our cosmology, 1 arcmin = 0.2964 Mpc, so this is 1.186 to 7.113 Mpc.\n", + "* **Model:** NFW profile with $M_{200c}$ (mass inside the radius where the mean density is 200 times\n", + " the critical density) and a fixed concentration $c_{200c} = 4$.\n", + "* **Cosmology:** flat $\\Lambda$CDM with $H_0 = 70$ km/s/Mpc and $\\Omega_m = 0.3$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "cluster_z = 0.35\n", + "\n", + "# Galaxy selection (same as tutorial 304.3)\n", + "band_list = [\"r\", \"i\", \"z\"]\n", + "snr_min = 5.0\n", + "color_rs = {\"ri\": 0.6, \"iz\": 0.3} # Red-sequence colors [AB mag]\n", + "rs_half_width = 0.1 # Half width of the red sequence [AB mag]\n", + "bg_mag_min = 18.0 # Background galaxies must be fainter than this [AB mag]\n", + "rs_mag_max = 20.0 # Red-sequence galaxies shown for reference are brighter than this [AB mag]\n", + "\n", + "# Shape-quality cuts (same as tutorial 304.3)\n", + "e_max = 2.0\n", + "sigma_e_max = 0.4\n", + "blendedness_max = 0.42\n", + "\n", + "# Shear calibration\n", + "e_rms = 0.365 # Assumed per-component intrinsic shape dispersion\n", + "m_bias = 0.0 # Multiplicative shear bias\n", + "c1_bias = 0.0 # Additive shear bias (component 1)\n", + "c2_bias = 0.0 # Additive shear bias (component 2)\n", + "\n", + "# Radial profile: logarithmic bins in projected physical distance [Mpc]\n", + "# Default: the tutorial 304.3 bins, 4 to 24 arcmin (~ 1 to 7 Mpc at z = 0.35) in 5 bins\n", + "r_min_mpc = 1.2\n", + "r_max_mpc = 7.1\n", + "n_bins = 5\n", + "\n", + "# Halo model\n", + "concentration = 4.0\n", + "massdef = \"critical\"\n", + "delta_mdef = 200\n", + "\n", + "# Assumed source redshift distribution (no photo-z)\n", + "z_distrib_func = redshift_distributions.chang2013 # Alternative: redshift_distributions.desc_srd\n", + "z_src_min = cluster_z + 0.1 # Sources are assumed to be behind the cluster\n", + "z_src_max = 5.0\n", + "z_inf = 10.0 # \"Infinite\" redshift used to normalize the lensing efficiency\n", + "\n", + "# Cosmology\n", + "cosmo = clmm.Cosmology(H0=70.0, Omega_dm0=0.3 - 0.045, Omega_b0=0.045, Omega_k0=0.0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 4. Select galaxies\n", + "\n", + "Keep objects that are extended (`refExtendedness == 1`, i.e. not point-like stars), have no\n", + "CModel failure flag, and have SNR > 5 in all of r, i and z." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sel = obj_cat[\"refExtendedness\"] == 1\n", + "for band in band_list:\n", + " sel &= obj_cat[f\"{band}_cModel_flag\"] == 0\n", + " sel &= obj_cat[f\"{band}_cModel_SNR\"] > snr_min\n", + "gal_cat = obj_cat[sel] # New table: obj_cat itself is never modified\n", + "print(\"Number of galaxies is\", len(gal_cat))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 5. Select background galaxies by color\n", + "\n", + "Cluster galaxies form a *red sequence* (RS) in the color–magnitude diagram (CMD).\n", + "Galaxies that are **redder than the RS** in both $r-i$ and $i-z$ are likely to be at higher redshift,\n", + "i.e. behind the cluster. We use them as lensing sources. Red-sequence galaxies (bright, inside the RS)\n", + "are selected only to show where the cluster is." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "color_ri = gal_cat[\"r_cModelMag\"] - gal_cat[\"i_cModelMag\"]\n", + "color_iz = gal_cat[\"i_cModelMag\"] - gal_cat[\"z_cModelMag\"]\n", + "\n", + "sel_bg = (color_ri > color_rs[\"ri\"] + rs_half_width) & (color_iz > color_rs[\"iz\"] + rs_half_width)\n", + "sel_rs = (np.abs(color_ri - color_rs[\"ri\"]) < rs_half_width) & (\n", + " np.abs(color_iz - color_rs[\"iz\"]) < rs_half_width\n", + ")\n", + "for band in band_list:\n", + " sel_bg &= gal_cat[f\"{band}_cModelMag\"] > bg_mag_min\n", + " sel_rs &= gal_cat[f\"{band}_cModelMag\"] < rs_mag_max\n", + "\n", + "print(\"Number of background galaxies is\", np.sum(sel_bg))\n", + "print(\"Number of red-sequence galaxies is\", np.sum(sel_rs))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot the CMDs. Dashed lines mark the RS; blue points are background galaxies and red points are\n", + "RS galaxies." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(1, 2, figsize=(12, 5))\n", + "for i, (band1, band2, color) in enumerate([(\"r\", \"i\", color_ri), (\"i\", \"z\", color_iz)]):\n", + " mag = gal_cat[f\"{band1}_cModelMag\"]\n", + " ax[i].scatter(mag, color, marker=\".\", s=0.3, alpha=0.5, color=\"k\", label=\"All galaxies\")\n", + " ax[i].scatter(mag[sel_bg], color[sel_bg], marker=\".\", s=0.3, color=\"blue\", label=\"Background\")\n", + " ax[i].scatter(mag[sel_rs], color[sel_rs], marker=\".\", s=0.3, color=\"red\", label=\"Red sequence\")\n", + " rs_color = color_rs[band1 + band2]\n", + " for offset in (-rs_half_width, rs_half_width):\n", + " ax[i].axhline(rs_color + offset, ls=\"--\", color=\"grey\")\n", + " ax[i].set_xlabel(f\"{band1} [AB mag]\")\n", + " ax[i].set_ylabel(f\"{band1}-{band2} [AB mag]\")\n", + " ax[i].set_xlim(17, 25)\n", + " ax[i].set_ylim(-0.5, 1.5)\n", + " ax[i].legend(markerscale=20, fontsize=0.8*label_fontsize)\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot the sky positions. RS galaxies (red) are concentrated around the cluster, while background\n", + "galaxies (blue) are spread nearly uniformly." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(6.5, 6.5))\n", + "ax.scatter(gal_cat[\"coord_ra\"][sel_bg], gal_cat[\"coord_dec\"][sel_bg], s=1, alpha=0.2, color=\"blue\")\n", + "ax.scatter(gal_cat[\"coord_ra\"][sel_rs], gal_cat[\"coord_dec\"][sel_rs], s=3, color=\"red\")\n", + "ax.invert_xaxis()\n", + "ax.set_xlabel(\"RA [deg]\")\n", + "ax.set_ylabel(\"Dec [deg]\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 6. Apply shape-quality cuts\n", + "\n", + "Following tutorial 304.3 (and Mandelbaum et al. 2018), keep background galaxies with:\n", + "* total ellipticity $\\sqrt{e_1^2 + e_2^2} < 2$ (removes unphysical shapes),\n", + "* shape uncertainty `hsmShapeRegauss_sigma` ≤ 0.4,\n", + "* no HSM failure flag,\n", + "* blendedness < 0.42 (flux not strongly affected by neighbors)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "e1_col = f\"{band_shape}_hsmShapeRegauss_e1\"\n", + "e2_col = f\"{band_shape}_hsmShapeRegauss_e2\"\n", + "\n", + "sel_shape = np.sqrt(gal_cat[e1_col] ** 2 + gal_cat[e2_col] ** 2) < e_max\n", + "sel_shape &= gal_cat[f\"{band_shape}_hsmShapeRegauss_sigma\"] <= sigma_e_max\n", + "sel_shape &= gal_cat[f\"{band_shape}_hsmShapeRegauss_flag\"] == 0\n", + "sel_shape &= gal_cat[f\"{band_shape}_blendedness\"] < blendedness_max\n", + "\n", + "src_cat = gal_cat[sel_bg & sel_shape]\n", + "print(\"Number of source galaxies is\", len(src_cat))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 7. From ellipticity to shear\n", + "\n", + "The HSM Regaussianization ellipticity is a *distortion*,\n", + "$e = (a^2 - b^2)/(a^2 + b^2)$ for an ellipse with axes $a$ and $b$. For a population of galaxies, the\n", + "mean distortion responds to a shear $g$ as $\\langle e \\rangle \\approx 2\\mathcal{R}\\, g$, where the\n", + "*responsivity* is\n", + "\n", + "$$\n", + "\\mathcal{R} = 1 - e_{\\rm rms}^2 ,\n", + "$$\n", + "\n", + "and $e_{\\rm rms}$ is the per-component intrinsic shape dispersion. Including multiplicative ($m$) and\n", + "additive ($c$) biases, the shear estimate of each galaxy is (Mandelbaum et al. 2018):\n", + "\n", + "$$\n", + "g_{1,2} = \\frac{1}{1+m}\\left(\\frac{e_{1,2}}{2\\mathcal{R}} - c_{1,2}\\right).\n", + "$$\n", + "\n", + "Here $e_{\\rm rms} = 0.365$ and $m = c = 0$ (see Section 3)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "responsivity = 1.0 - e_rms**2\n", + "print(f\"Responsivity R = {responsivity:.3f}\")\n", + "\n", + "e1 = np.asarray(src_cat[e1_col], dtype=float)\n", + "e2 = np.asarray(src_cat[e2_col], dtype=float)\n", + "g1 = (e1 / (2.0 * responsivity) - c1_bias) / (1.0 + m_bias)\n", + "g2 = (e2 / (2.0 * responsivity) - c2_bias) / (1.0 + m_bias)\n", + "\n", + "src_ra = np.asarray(src_cat[\"coord_ra\"], dtype=float)\n", + "src_dec = np.asarray(src_cat[\"coord_dec\"], dtype=float)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 8. Tangential and cross shear with CLMM\n", + "\n", + "Lensing stretches background galaxies *tangentially* around the cluster. For each source, CLMM\n", + "computes the angular separation $\\theta$ from the cluster center, the tangential shear\n", + "$g_t = -(g_1\\cos 2\\varphi + g_2 \\sin 2\\varphi)$ and the cross shear\n", + "$g_\\times = g_1 \\sin 2\\varphi - g_2\\cos 2\\varphi$, where $\\varphi$ is the position angle of the source.\n", + "\n", + "DP2 shapes use the same convention as CLMM's `coordinate_system=\"euclidean\"`\n", + "($\\varphi$ measured from West, counterclockwise), so no sign change is needed.\n", + "`geometry=\"curve\"` uses exact spherical geometry." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "theta, gt, gx = da.compute_tangential_and_cross_components(\n", + " ra_lens=cluster_ra,\n", + " dec_lens=cluster_dec,\n", + " ra_source=src_ra,\n", + " dec_source=src_dec,\n", + " shear1=g1,\n", + " shear2=g2,\n", + " coordinate_system=\"euclidean\",\n", + " geometry=\"curve\",\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Check:** compare with the formula used in tutorial 304.3 (astropy position angle + $\\pi/2$).\n", + "The two must agree." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "center = SkyCoord(cluster_ra, cluster_dec, frame=\"icrs\", unit=\"deg\")\n", + "sources = SkyCoord(src_ra, src_dec, frame=\"icrs\", unit=\"deg\")\n", + "phi = center.position_angle(sources).rad + np.pi / 2.0\n", + "gt_tutorial = -g1 * np.cos(2.0 * phi) - g2 * np.sin(2.0 * phi)\n", + "gx_tutorial = g1 * np.sin(2.0 * phi) - g2 * np.cos(2.0 * phi)\n", + "\n", + "print(\"Tangential shear agrees:\", np.allclose(gt, gt_tutorial))\n", + "print(\"Cross shear agrees:\", np.allclose(gx, gx_tutorial))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 9. Radial shear profile\n", + "\n", + "Average $g_t$ and $g_\\times$ in logarithmic bins of projected physical distance $R$ (Mpc), using the\n", + "cluster redshift and our cosmology to convert angles into distances. The error bar in each bin is the\n", + "standard error of the mean ($\\sigma/\\sqrt{N}$), dominated by galaxy shape noise.\n", + "\n", + "By default we use the bins of tutorial 304.3 (5 logarithmic bins from 4 to 24 arcmin, i.e. 1.186 to\n", + "7.113 Mpc). The center is skipped because there the source sample is most contaminated by cluster\n", + "members and the data are noisy. The next cell prints the bin edges in both Mpc and arcmin, so you\n", + "can check the conversion." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "bin_edges = make_bins(r_min_mpc, r_max_mpc, n_bins, method=\"evenlog10width\")\n", + "bin_edges_arcmin = clmm.utils.convert_units(bin_edges, \"Mpc\", \"arcmin\", cluster_z, cosmo)\n", + "print(\"Bin edges [Mpc]: \", np.round(bin_edges, 3))\n", + "print(\"Bin edges [arcmin]:\", np.round(bin_edges_arcmin, 2))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Compute the binned profile." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "profile = da.make_radial_profile(\n", + " [gt, gx],\n", + " angsep=theta,\n", + " angsep_units=\"radians\",\n", + " bin_units=\"Mpc\",\n", + " bins=bin_edges,\n", + " cosmo=cosmo,\n", + " z_lens=cluster_z,\n", + " coordinate_system=\"euclidean\",\n", + ")\n", + "# make_radial_profile names the components p_0, p_1; give them clearer names.\n", + "profile.rename_columns([\"p_0\", \"p_0_err\", \"p_1\", \"p_1_err\"], [\"gt\", \"gt_err\", \"gx\", \"gx_err\"])\n", + "profile[\"radius\", \"gt\", \"gt_err\", \"gx\", \"gx_err\", \"n_src\"]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot the profile. A lensing signal shows up as positive $g_t$ that decreases with radius, while\n", + "$g_\\times$ should be consistent with zero." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(8, 5.5))\n", + "ax.errorbar(profile[\"radius\"], profile[\"gt\"], profile[\"gt_err\"], fmt=\"o\", label=r\"$g_t$\")\n", + "ax.errorbar(\n", + " profile[\"radius\"] * 1.03, profile[\"gx\"], profile[\"gx_err\"], fmt=\"x\", label=r\"$g_\\times$\"\n", + ")\n", + "ax.axhline(0, color=\"k\", ls=\":\")\n", + "ax.set_xscale(\"log\")\n", + "ax.set_xlabel(\"R [Mpc]\")\n", + "ax.set_ylabel(\"Reduced shear\")\n", + "ax.legend(fontsize=0.8*label_fontsize)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 10. Assume a redshift distribution for the sources\n", + "\n", + "The lensing signal depends on how far the sources are behind the cluster, through the lensing\n", + "efficiency\n", + "\n", + "$$\n", + "\\beta_s(z_s) = \\frac{D_{LS}}{D_S}\\,\\frac{D_\\infty}{D_{L\\infty}},\n", + "$$\n", + "\n", + "where $D_{LS}$ and $D_S$ are angular-diameter distances from the lens and from us to the source.\n", + "Without photo-z, we assume a source redshift distribution $N(z)$, here the Chang et al. (2013)\n", + "distribution in CLMM, restricted to $z > z_{\\rm cluster} + 0.1$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "z_grid = np.linspace(0.0, z_src_max, 500)\n", + "norm = quad(z_distrib_func, 0.0, z_src_max)[0]\n", + "\n", + "fig, ax = plt.subplots(figsize=(8, 5))\n", + "ax.plot(z_grid, z_distrib_func(z_grid) / norm, label=\"Assumed N(z)\")\n", + "ax.axvline(cluster_z, color=\"orange\", label=\"Cluster redshift\")\n", + "ax.axvspan(0.0, z_src_min, color=\"grey\", alpha=0.2, label=\"Excluded (z < z_cl + 0.1)\")\n", + "ax.set_xlabel(\"Source redshift z\")\n", + "ax.set_ylabel(\"N(z) (normalized)\")\n", + "ax.legend(fontsize=0.8*label_fontsize)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Average $\\beta_s$ and $\\beta_s^2$ over $N(z)$. These two numbers are all the model needs to know about\n", + "the sources." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "beta_kwargs = {\n", + " \"z_cl\": cluster_z,\n", + " \"z_inf\": z_inf,\n", + " \"cosmo\": cosmo,\n", + " \"zmin\": z_src_min,\n", + " \"zmax\": z_src_max,\n", + " \"z_distrib_func\": z_distrib_func,\n", + "}\n", + "beta_s_mean = compute_beta_s_mean_from_distribution(**beta_kwargs)\n", + "beta_s_square_mean = compute_beta_s_square_mean_from_distribution(**beta_kwargs)\n", + "print(f\" = {beta_s_mean:.3f}\")\n", + "print(f\" = {beta_s_square_mean:.3f}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 11. Fit the halo mass\n", + "\n", + "The model is the reduced tangential shear of an NFW halo with mass $M_{200c}$ and fixed concentration.\n", + "With `z_src_info=\"beta\"` and `approx=\"type1\"`, CLMM uses the approximation of Applegate et al. (2014):\n", + "\n", + "$$\n", + "g_t(R) \\approx \\frac{\\langle\\beta_s\\rangle\\,\\gamma_\\infty(R)}\n", + "{1 - \\frac{\\langle\\beta_s^2\\rangle}{\\langle\\beta_s\\rangle}\\kappa_\\infty(R)} ,\n", + "$$\n", + "\n", + "where $\\gamma_\\infty$ and $\\kappa_\\infty$ are the shear and convergence for a source at $z_\\infty$.\n", + "We fit $\\log_{10} M_{200c}$ with `scipy.optimize.curve_fit` (through CLMM's `fitters`)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def gt_model(radius, logm):\n", + " return clmm.compute_reduced_tangential_shear(\n", + " r_proj=radius,\n", + " mdelta=10**logm,\n", + " cdelta=concentration,\n", + " z_cluster=cluster_z,\n", + " z_src=(beta_s_mean, beta_s_square_mean),\n", + " cosmo=cosmo,\n", + " delta_mdef=delta_mdef,\n", + " massdef=massdef,\n", + " halo_profile_model=\"nfw\",\n", + " z_src_info=\"beta\",\n", + " approx=\"type1\",\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "popt, pcov = fitters[\"curve_fit\"](\n", + " gt_model,\n", + " profile[\"radius\"],\n", + " profile[\"gt\"],\n", + " profile[\"gt_err\"],\n", + " bounds=[12.0, 16.0],\n", + ")\n", + "logm_fit, logm_err = popt[0], np.sqrt(pcov[0][0])\n", + "m_fit = 10**logm_fit\n", + "m_err = m_fit * np.log(10) * logm_err\n", + "\n", + "print(f\"log10(M200c / Msun) = {logm_fit:.2f} +/- {logm_err:.2f}\")\n", + "print(f\"M200c = ({m_fit / 1e14:.2f} +/- {m_err / 1e14:.2f}) x 10^14 Msun\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 12. Plot the result\n", + "\n", + "Data points with the best-fit NFW model. The shaded band shows the model for $\\log_{10} M_{200c} \\pm 1\\sigma$, where $1\\sigma$ is the statistical uncertainty of the fitted $\\log_{10} M_{200c}$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "r_model = np.logspace(np.log10(r_min_mpc), np.log10(r_max_mpc), 100)\n", + "\n", + "fig, ax = plt.subplots(figsize=(9, 6))\n", + "ax.errorbar(profile[\"radius\"], profile[\"gt\"], profile[\"gt_err\"], fmt=\"o\", c=\"k\", label=r\"$g_t$ data\")\n", + "ax.errorbar(\n", + " profile[\"radius\"] * 1.03,\n", + " profile[\"gx\"],\n", + " profile[\"gx_err\"],\n", + " fmt=\"x\",\n", + " c=\"grey\",\n", + " label=r\"$g_\\times$ data\",\n", + ")\n", + "ax.plot(\n", + " r_model,\n", + " gt_model(r_model, logm_fit),\n", + " \"C1\",\n", + " label=rf\"NFW fit: $M_{{200c}} = ({m_fit/1e14:.1f} \\pm {m_err/1e14:.1f})\\times10^{{14}}\\,M_\\odot$\",\n", + ")\n", + "ax.fill_between(\n", + " r_model,\n", + " gt_model(r_model, logm_fit - logm_err),\n", + " gt_model(r_model, logm_fit + logm_err),\n", + " color=\"C1\",\n", + " alpha=0.3,\n", + ")\n", + "ax.axhline(0, color=\"k\", ls=\":\")\n", + "ax.set_xscale(\"log\")\n", + "ax.set_xlabel(\"R [Mpc]\")\n", + "ax.set_ylabel(\"Reduced shear\")\n", + "ax.legend(loc=\"upper right\", fontsize=0.8*label_fontsize)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 13. Caveats and exercises\n", + "\n", + "This is a simple demonstration. The main simplifications are:\n", + "\n", + "* **Contamination by cluster members.** A color cut does not remove all cluster galaxies. Members\n", + " carry no lensing signal, so they dilute $g_t$ and bias the mass **low**, especially near the center.\n", + " A boost-factor correction can account for this (see `examples/demo_boost_factors.ipynb` in CLMM).\n", + "* **Assumed $N(z)$.** The redshift distribution of color-selected, relatively bright DP2 galaxies is\n", + " not measured here. A different $N(z)$ changes $\\langle\\beta_s\\rangle$ and therefore the mass.\n", + "* **Shear calibration.** We set $m = c = 0$ and fixed $e_{\\rm rms} = 0.365$; galaxies are not weighted.\n", + "* **Model.** Fixed concentration, BCG as the center (no miscentering), and no 2-halo term. The\n", + " default bins (tutorial 304.3) reach 7.1 Mpc, where the 2-halo term from surrounding large-scale\n", + " structure is no longer negligible. The field also contains another cluster to the south (see the\n", + " sky-position plot), which can affect the outer bins.\n", + "\n", + "**Exercises** (each needs only Section 3 and later to be rerun; Exercise 1 also needs Section 2)\n", + "\n", + "1. Use the r band for shapes (`band_shape = \"r\"` in Section 2) and compare the mass.\n", + "2. Use `redshift_distributions.desc_srd` for $N(z)$, or change `z_src_min`.\n", + "3. Use `approx=\"type2\"` (Schrabback et al. 2018) in `gt_model`.\n", + "4. Change the radial range (`r_min_mpc`, `r_max_mpc`) or `n_bins`, for example `r_max_mpc = 3.0`\n", + " to reduce the 2-halo contribution.\n", + "5. Replace the assumed $N(z)$ with photometric redshifts when they are available." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "LSST", + "language": "python", + "name": "lsst" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.9" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/sprint_week_2026/DP2_cluster_mass_fit_photoz-cleaned.ipynb b/sprint_week_2026/DP2_cluster_mass_fit_photoz-cleaned.ipynb new file mode 100644 index 0000000..d6085a1 --- /dev/null +++ b/sprint_week_2026/DP2_cluster_mass_fit_photoz-cleaned.ipynb @@ -0,0 +1,950 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Galaxy cluster weak-lensing mass with DP2 photo-z and CLMM\n", + "\n", + "For the Rubin Science Platform at data.lsst.cloud.
\n", + "Data Release: Data Preview 2
\n", + "Container Size: large
\n", + "LSST Science Pipelines version: r30 (same as DP2 tutorials 304.3 and 310.1)
\n", + "Packages: `lsst.rsp`, `lsdb`, `clmm` (≥ 1.18), a CLMM theory backend (`pyccl`)\n", + "\n", + "**Learning objective:** Select background galaxies behind the galaxy cluster PSZ2 G309.43-72.86 with\n", + "DP2 photometric redshifts (photo-z), measure their reduced tangential shear profile, and fit the\n", + "cluster mass with CLMM.\n", + "\n", + "**Credit:** Object catalog access and galaxy selection follow the Rubin DP2 tutorial\n", + "[304.3 \"Galaxy cluster weak lensing with HSM\"](https://dp2.lsst.io/tutorials/notebook/304/notebook-304-3.html);\n", + "photo-z access follows\n", + "[310.1 \"Photo-z access and visualization\"](https://dp2.lsst.io/tutorials/notebook/310/notebook-310-1.html).\n", + "This notebook extends `DP2_cluster_mass_fit_color_selection.ipynb` (in this folder), which selected\n", + "sources by color and assumed a redshift distribution." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Overview\n", + "\n", + "PSZ2 G309.43-72.86 is a rich galaxy cluster at redshift $z=0.35$ in DP2. In the previous notebook we\n", + "selected background galaxies by color and *assumed* their redshift distribution. Here we use the DP2\n", + "photo-z instead: a galaxy is a lensing source if its photo-z is larger than $z_{\\rm cluster} + 0.1$,\n", + "and the photo-z of the sources also give the lensing efficiency used in the mass model.\n", + "\n", + "1. **Setup** – make CLMM importable on the RSP.\n", + "2. **Load the data** (run only once): the DP2 Object catalog (TAP) and the photo-z (LSDB), matched by\n", + " `objectId`.\n", + "3. **Define the analysis parameters** in one place.\n", + "4. **Select galaxies** (extended, well-measured objects).\n", + "5. **Select background galaxies with photo-z.**\n", + "6. **Apply shape-quality cuts.**\n", + "7. **Convert HSM ellipticities into shear estimates.**\n", + "8. **Compute tangential and cross shear** with CLMM.\n", + "9. **Build the radial shear profile** in physical units (Mpc).\n", + "10. **Compute the lensing efficiency** from the source photo-z.\n", + "11. **Fit an NFW halo mass** with CLMM.\n", + "12. **Plot the result.**\n", + "13. **Caveats and exercises.**\n", + "\n", + "**Tip:** the catalogs from Section 2 are never modified afterwards. To try other choices (for example\n", + "another photo-z estimator, or the radial range), change the parameters in Section 3 and rerun from\n", + "Section 3 onwards, without loading the data again." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1. Setup\n", + "\n", + "### 1.1 One-time installation (only the first time)\n", + "\n", + "CLMM is not part of the LSST Science Pipelines, so we use the latest CLMM code from GitHub.\n", + "Open a terminal on the RSP, go to any directory you like, and clone CLMM:\n", + "\n", + "```bash\n", + "git clone https://github.com/LSSTDESC/CLMM.git\n", + "```\n", + "\n", + "To get newer CLMM code later, run `git pull` inside the `CLMM` folder.\n", + "\n", + "`import clmm` also needs `qp` (used by CLMM for photo-z PDFs; imported even though this notebook\n", + "does not use photo-z) and a theory backend (we use CCL). If they are not available in your\n", + "RSP environment, uncomment and run the cell below **once**, then restart the kernel\n", + "(Kernel → Restart Kernel)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# %pip install --user qp-prob pyccl" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 1.2 Make CLMM importable\n", + "\n", + "Jupyter runs a notebook from the notebook's own folder, and the CLMM clone's folder needs to be provided. We put it at the front of the Python path, so `import clmm` uses this\n", + "GitHub version instead of any older installed version. We also ask CLMM to use the CCL backend." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import sys\n", + "\n", + "CLMM_PATH = os.path.abspath(\"./CLMM\") # The clone's top folder, which contains the \"clmm\" package\n", + "sys.path.insert(0, CLMM_PATH)\n", + "os.environ[\"CLMM_MODELING_BACKEND\"] = \"ccl\" # Must be set before importing clmm\n", + "\n", + "import clmm\n", + "\n", + "print(\"CLMM version:\", clmm.__version__)\n", + "print(\"CLMM loaded from:\", clmm.__file__)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 1.3 Import packages\n", + "\n", + "General packages (numpy, scipy, matplotlib, astropy), the RSP data-access packages (`lsst.rsp` for\n", + "TAP and `lsdb` for the photo-z) and the CLMM modules we need." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import lsdb\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from astropy.coordinates import SkyCoord\n", + "from astropy.table import Table, join\n", + "from lsst.rsp import RSPDiscovery\n", + "from scipy.integrate import quad\n", + "\n", + "import clmm.dataops as da\n", + "from clmm.support.sampler import fitters\n", + "from clmm.utils import (\n", + " compute_beta_s,\n", + " compute_beta_s_mean_from_distribution,\n", + " compute_beta_s_square_mean_from_distribution,\n", + " make_bins,\n", + " redshift_distributions,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Make the axis labels and tick labels in all figures 1.7 times larger than the matplotlib default." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "label_fontsize = 1.7 * plt.rcParams[\"font.size\"]\n", + "plt.rcParams.update(\n", + " {\n", + " \"axes.labelsize\": label_fontsize,\n", + " \"xtick.labelsize\": label_fontsize,\n", + " \"ytick.labelsize\": label_fontsize,\n", + " }\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2. Load the data\n", + "\n", + "This section is run **once**. It depends only on the cluster position (we use the Brightest Cluster\n", + "Galaxy, BCG, as the center), the band used for galaxy shapes, and the search radius." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Cluster PSZ2 G309.43-72.86 (BCG position used as the center)\n", + "cluster_ra = 10.2082 # [deg]\n", + "cluster_dec = -44.1307 # [deg]\n", + "\n", + "band_shape = \"i\" # Band used for galaxy shapes\n", + "radius_deg = 0.4 # Cone-search radius around the BCG [deg]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 2.1 Object catalog (TAP)\n", + "\n", + "Get the TAP service, exactly as in tutorial 304.3." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "discovery = RSPDiscovery(\"dp2\")\n", + "service = discovery.get_tap_client()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Query the `dp2.Object` table for objects within `radius_deg` of the BCG. We ask for the HSM\n", + "Regaussianization shapes in the `band_shape` band, the riz CModel photometry, and quality flags.\n", + "The signal-to-noise ratio (SNR) is computed in the query because it is not a column in the schema." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "query = (\n", + " \"SELECT objectId, coord_ra, coord_dec, refExtendedness, \"\n", + " f\"{band_shape}_blendedness, \"\n", + " f\"{band_shape}_hsmShapeRegauss_e1, \"\n", + " f\"{band_shape}_hsmShapeRegauss_e2, \"\n", + " f\"{band_shape}_hsmShapeRegauss_sigma, \"\n", + " f\"{band_shape}_hsmShapeRegauss_flag, \"\n", + " \"r_cModelMag, i_cModelMag, z_cModelMag, \"\n", + " \"r_cModelFlux/r_cModelFluxErr AS r_cModel_SNR, \"\n", + " \"i_cModelFlux/i_cModelFluxErr AS i_cModel_SNR, \"\n", + " \"z_cModelFlux/z_cModelFluxErr AS z_cModel_SNR, \"\n", + " \"r_cModel_flag, i_cModel_flag, z_cModel_flag \"\n", + " \"FROM dp2.Object \"\n", + " \"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), \"\n", + " f\"CIRCLE('ICRS', {cluster_ra}, {cluster_dec}, {radius_deg})) = 1 \"\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Run the query, check that it completed, and convert the result into an astropy table." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "job = service.submit_job(query)\n", + "job.run()\n", + "job.wait(phases=[\"COMPLETED\", \"ERROR\"])\n", + "print(\"Job phase is\", job.phase)\n", + "if job.phase == \"ERROR\":\n", + " job.raise_if_error()\n", + "assert job.phase == \"COMPLETED\"\n", + "obj_cat = job.fetch_result().to_table()\n", + "print(\"Number of objects is\", len(obj_cat))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 2.2 Photo-z catalog (LSDB)\n", + "\n", + "The DP2 photo-z are provided as LSDB-formatted files on the RSP (not as a TAP table), described in\n", + "\"Photometric Redshifts for Data Preview 2\" (RTN-124). Six estimators are available: `bpz`, `dnf`,\n", + "`fzboost`, `gpz`, `knn` and `tpz`. For each, we load the point estimate `z_best` (the minimal-risk\n", + "estimate of Tanaka et al. 2018) and the lower and upper ends of the 68% credible interval.\n", + "\n", + "We load all six estimators at once (the cone is small), so that changing the estimator in Section 3\n", + "does not require loading the data again. Only the objects within `radius_deg` of the BCG are read.\n", + "\n", + "Note: the DP2 photo-z are only lightly validated and have no quality flags yet; treat them as\n", + "provisional (see tutorial 310.2)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "pz_codes = [\"bpz\", \"dnf\", \"fzboost\", \"gpz\", \"knn\", \"tpz\"]\n", + "pz_columns = []\n", + "for pz_code_i in pz_codes:\n", + " pz_columns += [\n", + " f\"{pz_code_i}_z_best\",\n", + " f\"{pz_code_i}_z_err68_low\",\n", + " f\"{pz_code_i}_z_err68_high\",\n", + " ]\n", + "\n", + "pz_lsdb = lsdb.open_catalog(\n", + " \"/rubin/lsdb_data/dp2/object_photoz\",\n", + " search_filter=lsdb.ConeSearch(\n", + " ra=cluster_ra, dec=cluster_dec, radius_arcsec=radius_deg * 3600.0\n", + " ),\n", + " columns=[\"objectId\", *pz_columns],\n", + ")\n", + "pz_df = pz_lsdb.compute()\n", + "print(\"Number of objects with photo-z is\", len(pz_df))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 2.3 Match the photo-z to the Object catalog\n", + "\n", + "Both catalogs share the `objectId` column, so we match them by ID (no positional cross-match is\n", + "needed). First convert the photo-z table into an astropy table with plain numpy columns (missing\n", + "values become NaN), then keep the objects present in both catalogs (an \"inner\" join)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "pz_table = Table()\n", + "pz_table[\"objectId\"] = pz_df[\"objectId\"].to_numpy(dtype=np.int64)\n", + "for col in pz_columns:\n", + " pz_table[col] = pz_df[col].to_numpy(dtype=float, na_value=np.nan)\n", + "\n", + "obj_pz_cat = join(obj_cat, pz_table, keys=\"objectId\", join_type=\"inner\")\n", + "print(\"Number of objects in the Object catalog:\", len(obj_cat))\n", + "print(\"Number of objects with a photo-z match: \", len(obj_pz_cat))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 3. Define the analysis parameters\n", + "\n", + "All analysis choices are collected here. **After changing any of them, rerun from this cell\n", + "onwards**; there is no need to load the data again.\n", + "\n", + "* **Cluster redshift:** $z = 0.35$.\n", + "* **Photo-z:** the FlexZBoost estimator (`fzboost`), which RTN-124 found to perform best on DP2 (see\n", + " tutorial 310.2). Sources must have photo-z $> z_{\\rm cluster} + 0.1$.\n", + "* **Selection:** the same galaxy and shape cuts as tutorial 304.3 (no red-sequence color cut).\n", + "* **Shear calibration:** DP2 does not provide multiplicative ($m$) or additive ($c$) shear biases\n", + " for the HSM shapes, so we set $m = c = 0$. We assume a per-component shape dispersion\n", + " $e_{\\rm rms} = 0.365$.\n", + "* **Radial bins:** the same bins as tutorial 304.3, which uses 5 logarithmic bins from 4 to 24 arcmin.\n", + " At $z = 0.35$ in our cosmology, 1 arcmin = 0.2964 Mpc, so this is 1.186 to 7.113 Mpc.\n", + "* **Model:** NFW profile with $M_{200c}$ (mass inside the radius where the mean density is 200 times\n", + " the critical density) and a fixed concentration $c_{200c} = 4$.\n", + "* **Cosmology:** flat $\\Lambda$CDM with $H_0 = 70$ km/s/Mpc and $\\Omega_m = 0.3$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "cluster_z = 0.35\n", + "\n", + "# Photo-z\n", + "pz_code = \"fzboost\" # Photo-z estimator: one of pz_codes (Section 2.2)\n", + "dz_min = 0.1 # Sources must have photo-z > cluster_z + dz_min\n", + "\n", + "# Galaxy selection (same as tutorial 304.3)\n", + "band_list = [\"r\", \"i\", \"z\"]\n", + "snr_min = 5.0\n", + "\n", + "# Shape-quality cuts (same as tutorial 304.3)\n", + "e_max = 2.0\n", + "sigma_e_max = 0.4\n", + "blendedness_max = 0.42\n", + "\n", + "# Shear calibration\n", + "e_rms = 0.365 # Assumed per-component intrinsic shape dispersion\n", + "m_bias = 0.0 # Multiplicative shear bias\n", + "c1_bias = 0.0 # Additive shear bias (component 1)\n", + "c2_bias = 0.0 # Additive shear bias (component 2)\n", + "\n", + "# Radial profile: logarithmic bins in projected physical distance [Mpc]\n", + "# Default: the tutorial 304.3 bins, 4 to 24 arcmin (= 1.186 to 7.113 Mpc at z = 0.35) in 5 bins\n", + "r_min_mpc = 1.2\n", + "r_max_mpc = 7.1\n", + "n_bins = 5\n", + "\n", + "# Halo model\n", + "concentration = 4.0\n", + "massdef = \"critical\"\n", + "delta_mdef = 200\n", + "\n", + "# Lensing efficiency\n", + "z_inf = 10.0 # \"Infinite\" redshift used to normalize the lensing efficiency\n", + "\n", + "# Cosmology\n", + "cosmo = clmm.Cosmology(H0=70.0, Omega_dm0=0.3 - 0.045, Omega_b0=0.045, Omega_k0=0.0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 4. Select galaxies\n", + "\n", + "Keep objects that are extended (`refExtendedness == 1`, i.e. not point-like stars), have no\n", + "CModel failure flag, and have SNR > 5 in all of r, i and z." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sel = obj_pz_cat[\"refExtendedness\"] == 1\n", + "for band in band_list:\n", + " sel &= obj_pz_cat[f\"{band}_cModel_flag\"] == 0\n", + " sel &= obj_pz_cat[f\"{band}_cModel_SNR\"] > snr_min\n", + "gal_cat = obj_pz_cat[sel] # New table: obj_pz_cat itself is never modified\n", + "print(\"Number of galaxies is\", len(gal_cat))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 5. Select background galaxies with photo-z\n", + "\n", + "A galaxy can only be lensed by the cluster if it is behind the cluster. We keep galaxies whose photo-z\n", + "point estimate is larger than $z_{\\rm cluster} + 0.1$. The margin of 0.1 reduces contamination by\n", + "cluster members and foreground galaxies whose photo-z scatter slightly above the cluster redshift.\n", + "Galaxies without a valid photo-z are not selected." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "z_col = f\"{pz_code}_z_best\"\n", + "z_phot = np.asarray(gal_cat[z_col], dtype=float)\n", + "z_src_min = cluster_z + dz_min\n", + "\n", + "sel_bg = z_phot > z_src_min\n", + "print(f\"Photo-z estimator: {pz_code}\")\n", + "print(\"Number of galaxies with a valid photo-z is\", np.sum(np.isfinite(z_phot)))\n", + "print(f\"Number of background galaxies (z > {z_src_min:.2f}) is\", np.sum(sel_bg))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot the photo-z distribution of all galaxies and of the selected background galaxies. The orange line\n", + "marks the cluster redshift and the dashed line the selection cut." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "z_bins = np.linspace(0.0, 3.0, 61)\n", + "\n", + "fig, ax = plt.subplots(figsize=(8, 5))\n", + "ax.hist(z_phot[np.isfinite(z_phot)], bins=z_bins, histtype=\"step\", color=\"k\", label=\"All galaxies\")\n", + "ax.hist(z_phot[sel_bg], bins=z_bins, histtype=\"stepfilled\", alpha=0.4, label=\"Background galaxies\")\n", + "ax.axvline(cluster_z, color=\"orange\", label=\"Cluster redshift\")\n", + "ax.axvline(z_src_min, color=\"grey\", ls=\"--\", label=r\"$z_{\\rm cluster} + 0.1$\")\n", + "ax.set_xlabel(f\"Photo-z ({pz_code} z_best)\")\n", + "ax.set_ylabel(\"Number of galaxies\")\n", + "ax.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 6. Apply shape-quality cuts\n", + "\n", + "Following tutorial 304.3 (and Mandelbaum et al. 2018), keep background galaxies with:\n", + "* total ellipticity $\\sqrt{e_1^2 + e_2^2} < 2$ (removes unphysical shapes),\n", + "* shape uncertainty `hsmShapeRegauss_sigma` ≤ 0.4,\n", + "* no HSM failure flag,\n", + "* blendedness < 0.42 (flux not strongly affected by neighbors)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "e1_col = f\"{band_shape}_hsmShapeRegauss_e1\"\n", + "e2_col = f\"{band_shape}_hsmShapeRegauss_e2\"\n", + "\n", + "sel_shape = np.sqrt(gal_cat[e1_col] ** 2 + gal_cat[e2_col] ** 2) < e_max\n", + "sel_shape &= gal_cat[f\"{band_shape}_hsmShapeRegauss_sigma\"] <= sigma_e_max\n", + "sel_shape &= gal_cat[f\"{band_shape}_hsmShapeRegauss_flag\"] == 0\n", + "sel_shape &= gal_cat[f\"{band_shape}_blendedness\"] < blendedness_max\n", + "\n", + "src_cat = gal_cat[sel_bg & sel_shape]\n", + "print(\"Number of source galaxies is\", len(src_cat))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 7. From ellipticity to shear\n", + "\n", + "The HSM Regaussianization ellipticity is a *distortion*,\n", + "$e = (a^2 - b^2)/(a^2 + b^2)$ for an ellipse with axes $a$ and $b$. For a population of galaxies, the\n", + "mean distortion responds to a shear $g$ as $\\langle e \\rangle \\approx 2\\mathcal{R}\\, g$, where the\n", + "*responsivity* is\n", + "\n", + "$$\n", + "\\mathcal{R} = 1 - e_{\\rm rms}^2 ,\n", + "$$\n", + "\n", + "and $e_{\\rm rms}$ is the per-component intrinsic shape dispersion. Including multiplicative ($m$) and\n", + "additive ($c$) biases, the shear estimate of each galaxy is (Mandelbaum et al. 2018):\n", + "\n", + "$$\n", + "g_{1,2} = \\frac{1}{1+m}\\left(\\frac{e_{1,2}}{2\\mathcal{R}} - c_{1,2}\\right).\n", + "$$\n", + "\n", + "Here $e_{\\rm rms} = 0.365$ and $m = c = 0$ (see Section 3)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "responsivity = 1.0 - e_rms**2\n", + "print(f\"Responsivity R = {responsivity:.3f}\")\n", + "\n", + "e1 = np.asarray(src_cat[e1_col], dtype=float)\n", + "e2 = np.asarray(src_cat[e2_col], dtype=float)\n", + "g1 = (e1 / (2.0 * responsivity) - c1_bias) / (1.0 + m_bias)\n", + "g2 = (e2 / (2.0 * responsivity) - c2_bias) / (1.0 + m_bias)\n", + "\n", + "src_ra = np.asarray(src_cat[\"coord_ra\"], dtype=float)\n", + "src_dec = np.asarray(src_cat[\"coord_dec\"], dtype=float)\n", + "src_z = np.asarray(src_cat[z_col], dtype=float)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 8. Tangential and cross shear with CLMM\n", + "\n", + "Lensing stretches background galaxies *tangentially* around the cluster. For each source, CLMM\n", + "computes the angular separation $\\theta$ from the cluster center, the tangential shear\n", + "$g_t = -(g_1\\cos 2\\varphi + g_2 \\sin 2\\varphi)$ and the cross shear\n", + "$g_\\times = g_1 \\sin 2\\varphi - g_2\\cos 2\\varphi$, where $\\varphi$ is the position angle of the source.\n", + "\n", + "DP2 shapes use the same convention as CLMM's `coordinate_system=\"euclidean\"`\n", + "($\\varphi$ measured from West, counterclockwise), so no sign change is needed.\n", + "`geometry=\"curve\"` uses exact spherical geometry." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "theta, gt, gx = da.compute_tangential_and_cross_components(\n", + " ra_lens=cluster_ra,\n", + " dec_lens=cluster_dec,\n", + " ra_source=src_ra,\n", + " dec_source=src_dec,\n", + " shear1=g1,\n", + " shear2=g2,\n", + " coordinate_system=\"euclidean\",\n", + " geometry=\"curve\",\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Check:** compare with the formula used in tutorial 304.3 (astropy position angle + $\\pi/2$).\n", + "The two must agree." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "center = SkyCoord(cluster_ra, cluster_dec, frame=\"icrs\", unit=\"deg\")\n", + "sources = SkyCoord(src_ra, src_dec, frame=\"icrs\", unit=\"deg\")\n", + "phi = center.position_angle(sources).rad + np.pi / 2.0\n", + "gt_tutorial = -g1 * np.cos(2.0 * phi) - g2 * np.sin(2.0 * phi)\n", + "gx_tutorial = g1 * np.sin(2.0 * phi) - g2 * np.cos(2.0 * phi)\n", + "\n", + "print(\"Tangential shear agrees:\", np.allclose(gt, gt_tutorial))\n", + "print(\"Cross shear agrees:\", np.allclose(gx, gx_tutorial))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 9. Radial shear profile\n", + "\n", + "Average $g_t$ and $g_\\times$ in logarithmic bins of projected physical distance $R$ (Mpc), using the\n", + "cluster redshift and our cosmology to convert angles into distances. The error bar in each bin is the\n", + "standard error of the mean ($\\sigma/\\sqrt{N}$), dominated by galaxy shape noise.\n", + "\n", + "By default we use the bins of tutorial 304.3 (5 logarithmic bins from 4 to 24 arcmin, i.e. 1.186 to\n", + "7.113 Mpc). The center is skipped because there the source sample is most contaminated by cluster\n", + "members and the data are noisy. The next cell prints the bin edges in both Mpc and arcmin, so you\n", + "can check the conversion." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "bin_edges = make_bins(r_min_mpc, r_max_mpc, n_bins, method=\"evenlog10width\")\n", + "bin_edges_arcmin = clmm.utils.convert_units(bin_edges, \"Mpc\", \"arcmin\", cluster_z, cosmo)\n", + "print(\"Bin edges [Mpc]: \", np.round(bin_edges, 3))\n", + "print(\"Bin edges [arcmin]:\", np.round(bin_edges_arcmin, 2))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Compute the binned profile." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "profile = da.make_radial_profile(\n", + " [gt, gx],\n", + " angsep=theta,\n", + " angsep_units=\"radians\",\n", + " bin_units=\"Mpc\",\n", + " bins=bin_edges,\n", + " cosmo=cosmo,\n", + " z_lens=cluster_z,\n", + " coordinate_system=\"euclidean\",\n", + ")\n", + "# make_radial_profile names the components p_0, p_1; give them clearer names.\n", + "profile.rename_columns([\"p_0\", \"p_0_err\", \"p_1\", \"p_1_err\"], [\"gt\", \"gt_err\", \"gx\", \"gx_err\"])\n", + "profile[\"radius\", \"gt\", \"gt_err\", \"gx\", \"gx_err\", \"n_src\"]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot the profile. A lensing signal shows up as positive $g_t$ that decreases with radius, while\n", + "$g_\\times$ should be consistent with zero." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(8, 5.5))\n", + "ax.errorbar(profile[\"radius\"], profile[\"gt\"], profile[\"gt_err\"], fmt=\"o\", label=r\"$g_t$\")\n", + "ax.errorbar(\n", + " profile[\"radius\"] * 1.03, profile[\"gx\"], profile[\"gx_err\"], fmt=\"x\", label=r\"$g_\\times$\"\n", + ")\n", + "ax.axhline(0, color=\"k\", ls=\":\")\n", + "ax.set_xscale(\"log\")\n", + "ax.set_xlabel(\"R [Mpc]\")\n", + "ax.set_ylabel(\"Reduced shear\")\n", + "ax.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 10. Compute the lensing efficiency from the photo-z\n", + "\n", + "The lensing signal depends on how far each source is behind the cluster, through the lensing\n", + "efficiency\n", + "\n", + "$$\n", + "\\beta_s(z_s) = \\frac{D_{LS}}{D_S}\\,\\frac{D_\\infty}{D_{L\\infty}},\n", + "$$\n", + "\n", + "where $D_{LS}$ and $D_S$ are angular-diameter distances from the lens and from us to the source.\n", + "With photo-z, we compute $\\beta_s$ for every source and average $\\beta_s$ and $\\beta_s^2$ over the\n", + "source sample. These two numbers are all the model needs to know about the sources." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "beta_s = compute_beta_s(src_z, cluster_z, z_inf, cosmo)\n", + "beta_s_mean = np.mean(beta_s)\n", + "beta_s_square_mean = np.mean(beta_s**2)\n", + "print(f\" = {beta_s_mean:.3f}\")\n", + "print(f\" = {beta_s_square_mean:.3f}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For comparison, the previous notebook assumed the Chang et al. (2013) redshift distribution for the\n", + "sources. Compute its $\\langle\\beta_s\\rangle$ and $\\langle\\beta_s^2\\rangle$ with the same lower redshift\n", + "cut, and compare the two distributions." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "chang_kwargs = {\n", + " \"z_cl\": cluster_z,\n", + " \"z_inf\": z_inf,\n", + " \"cosmo\": cosmo,\n", + " \"zmin\": z_src_min,\n", + " \"zmax\": 5.0,\n", + " \"z_distrib_func\": redshift_distributions.chang2013,\n", + "}\n", + "print(f\"Chang+2013: = {compute_beta_s_mean_from_distribution(**chang_kwargs):.3f}\")\n", + "print(f\"Chang+2013: = {compute_beta_s_square_mean_from_distribution(**chang_kwargs):.3f}\")\n", + "\n", + "z_grid = np.linspace(z_src_min, 3.0, 300)\n", + "chang_norm = quad(redshift_distributions.chang2013, z_src_min, 5.0)[0]\n", + "\n", + "fig, ax = plt.subplots(figsize=(8, 5))\n", + "ax.hist(src_z, bins=z_bins, density=True, histtype=\"step\", color=\"k\", label=\"Sources (photo-z)\")\n", + "ax.plot(z_grid, redshift_distributions.chang2013(z_grid) / chang_norm, label=\"Chang+2013\")\n", + "ax.axvline(cluster_z, color=\"orange\", label=\"Cluster redshift\")\n", + "ax.set_xlabel(\"Source redshift z\")\n", + "ax.set_ylabel(\"N(z) (normalized)\")\n", + "ax.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 11. Fit the halo mass\n", + "\n", + "The model is the reduced tangential shear of an NFW halo with mass $M_{200c}$ and fixed concentration.\n", + "With `z_src_info=\"beta\"` and `approx=\"type1\"`, CLMM uses the approximation of Applegate et al. (2014):\n", + "\n", + "$$\n", + "g_t(R) \\approx \\frac{\\langle\\beta_s\\rangle\\,\\gamma_\\infty(R)}\n", + "{1 - \\frac{\\langle\\beta_s^2\\rangle}{\\langle\\beta_s\\rangle}\\kappa_\\infty(R)} ,\n", + "$$\n", + "\n", + "where $\\gamma_\\infty$ and $\\kappa_\\infty$ are the shear and convergence for a source at $z_\\infty$.\n", + "Here $\\langle\\beta_s\\rangle$ and $\\langle\\beta_s^2\\rangle$ come from the source photo-z (Section 10).\n", + "We fit $\\log_{10} M_{200c}$ with `scipy.optimize.curve_fit` (through CLMM's `fitters`)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def gt_model(radius, logm):\n", + " return clmm.compute_reduced_tangential_shear(\n", + " r_proj=radius,\n", + " mdelta=10**logm,\n", + " cdelta=concentration,\n", + " z_cluster=cluster_z,\n", + " z_src=(beta_s_mean, beta_s_square_mean),\n", + " cosmo=cosmo,\n", + " delta_mdef=delta_mdef,\n", + " massdef=massdef,\n", + " halo_profile_model=\"nfw\",\n", + " z_src_info=\"beta\",\n", + " approx=\"type1\",\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "popt, pcov = fitters[\"curve_fit\"](\n", + " gt_model,\n", + " profile[\"radius\"],\n", + " profile[\"gt\"],\n", + " profile[\"gt_err\"],\n", + " bounds=[12.0, 16.0],\n", + ")\n", + "logm_fit, logm_err = popt[0], np.sqrt(pcov[0][0])\n", + "m_fit = 10**logm_fit\n", + "m_err = m_fit * np.log(10) * logm_err\n", + "\n", + "print(f\"log10(M200c / Msun) = {logm_fit:.2f} +/- {logm_err:.2f}\")\n", + "print(f\"M200c = ({m_fit / 1e14:.2f} +/- {m_err / 1e14:.2f}) x 10^14 Msun\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 12. Plot the result\n", + "\n", + "Data points with the best-fit NFW model (orange line). The fit is done in $\\log_{10} M_{200c}$, and\n", + "$\\sigma$ is the 1$\\sigma$ statistical uncertainty of the fitted $\\log_{10} M_{200c}$ (from the\n", + "`curve_fit` covariance). The shaded band shows the model for $\\log_{10} M_{200c} \\pm \\sigma$, i.e. for\n", + "masses between $M_{200c} \\times 10^{-\\sigma}$ and $M_{200c} \\times 10^{+\\sigma}$ (concentration fixed).\n", + "\n", + "Note: the $\\pm$ on $M_{200c}$ in the legend is the first-order conversion\n", + "$\\sigma_M = M_{200c} \\ln(10)\\, \\sigma$, so it is symmetric, while the band is not." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "r_model = np.logspace(np.log10(r_min_mpc), np.log10(r_max_mpc), 100)\n", + "\n", + "fig, ax = plt.subplots(figsize=(9, 6))\n", + "ax.errorbar(profile[\"radius\"], profile[\"gt\"], profile[\"gt_err\"], fmt=\"o\", c=\"k\", label=r\"$g_t$ data\")\n", + "ax.errorbar(\n", + " profile[\"radius\"] * 1.03,\n", + " profile[\"gx\"],\n", + " profile[\"gx_err\"],\n", + " fmt=\"x\",\n", + " c=\"grey\",\n", + " label=r\"$g_\\times$ data\",\n", + ")\n", + "ax.plot(\n", + " r_model,\n", + " gt_model(r_model, logm_fit),\n", + " \"C1\",\n", + " label=rf\"NFW fit: $M_{{200c}} = ({m_fit/1e14:.1f} \\pm {m_err/1e14:.1f})\\times10^{{14}}\\,M_\\odot$\",\n", + ")\n", + "ax.fill_between(\n", + " r_model,\n", + " gt_model(r_model, logm_fit - logm_err),\n", + " gt_model(r_model, logm_fit + logm_err),\n", + " color=\"C1\",\n", + " alpha=0.3,\n", + ")\n", + "ax.axhline(0, color=\"k\", ls=\":\")\n", + "ax.set_xscale(\"log\")\n", + "ax.set_xlabel(\"R [Mpc]\")\n", + "ax.set_ylabel(\"Reduced shear\")\n", + "ax.legend(loc=\"upper right\", fontsize=label_fontsize)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 13. Caveats and exercises\n", + "\n", + "This is a simple demonstration. The main simplifications are:\n", + "\n", + "* **Photo-z errors.** Photo-z scatter and outliers put some cluster members and foreground galaxies\n", + " above the redshift cut. They carry no lensing signal (or a weaker one), so they dilute $g_t$ and\n", + " bias the mass **low**. Errors in the source photo-z also bias $\\langle\\beta_s\\rangle$, and hence the\n", + " mass. The DP2 photo-z are provisional and only lightly validated.\n", + "* **Point estimates.** We use one point estimate (`z_best`) per galaxy, not the full photo-z\n", + " probability distribution.\n", + "* **Shear calibration.** We set $m = c = 0$ and fixed $e_{\\rm rms} = 0.365$; galaxies are not weighted.\n", + "* **Model.** Fixed concentration, BCG as the center (no miscentering), and no 2-halo term. The\n", + " default bins (tutorial 304.3) reach 7.1 Mpc, where the 2-halo term from surrounding large-scale\n", + " structure is no longer negligible.\n", + "\n", + "**Exercises** (each needs only Section 3 and later to be rerun)\n", + "\n", + "1. Use the BPZ photo-z (`pz_code = \"bpz\"`) and compare the source sample and the mass.\n", + "2. Use a stricter source cut that also requires the lower end of the 68% photo-z interval to be\n", + " behind the cluster: `sel_bg = (z_phot > z_src_min) & (gal_cat[f\"{pz_code}_z_err68_low\"] > z_src_min)`.\n", + "3. Change `dz_min`, or the radial range (`r_min_mpc`, `r_max_mpc`).\n", + "4. Instead of $\\langle\\beta_s\\rangle$ and $\\langle\\beta_s^2\\rangle$, model each radial bin with the\n", + " individual source redshifts in that bin (`z_src_info=\"discrete\"`; see the \"z,R\" model in the CLMM\n", + " HSC example).\n", + "5. Compare the mass with the one from `DP2_cluster_mass_fit_color_selection.ipynb`." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "LSST", + "language": "python", + "name": "lsst" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.9" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/sprint_week_2026/readme.md b/sprint_week_2026/readme.md new file mode 100644 index 0000000..6539f00 --- /dev/null +++ b/sprint_week_2026/readme.md @@ -0,0 +1,44 @@ +# Readme + +We store the scripts here for the Sprint Week 2026. + +## [Oct 26, 2026] + +### Installing CLMM on RSP + +Run the code below: + +``` +pip install clmm +pip install pyccl +pip install qp-prob +``` + +Using `conda` may work, but it is under the lsst stack and may cause some issue. + +You can test the installation by running `import clmm` via `python3` in the terminal or in a Jupyter notebook. The current version is 1.16.10. + +To use the latest version of the code, use +``` +git clone https://github.com/LSSTDESC/CLMM.git +``` + +Then +``` +import sys +sys.path.insert(0, "path_to_CLMM") + +import clmm +print(clmm.__version__, clmm.__file__) # should print 1.18.0 (or the latest one) and a path inside your clone +``` + +Or equivalently +``` +cd path_to_CLMM +pip uninstall clmm # remove the 1.16.10 copy, if present +pip install -e . # links to your clone and installs the dependencies +pip install pyccl # or another backend +pip install qp-prob +``` + +This will use the current branch. Run `git pull` to get the newest `main`.