{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Day 3 - Get the Weak Lensing Signals" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this session, we will use the things described in the last two sessions and make a weak lensing measurement pipeline for ourselves. We will start with writing out the different steps required to get the pipeline ready. We then code up each of them as individual functions, which will be called in our main code.\n", "\n", "We are here studying the matter distribution around SDSS redmapper galaxy clusters - [arXiv:1303.3562](https://arxiv.org/abs/1303.3562) by the weak lensing technique using shape catalog data from HSC S16A data - [arXiv:1705.06745](https://arxiv.org/abs/1705.06745). We suggest readers to have a look at these references to have a better understanding of data provided to them. The readers can also compare their results with the same study done on these clusters using shape catalog from SDSS - [arXiv:1707.01907](https://arxiv.org/abs/1707.01907).\n", "\n", "**We highly suggest readers to make a separate notebook for this session as it will be useful for later use.** \n", "\n", "\n", "**In the codes below I have left some black spaces indicated using ??. Please fill up the ?? in the code below and make sure you are doing it with right equations and conversions.** " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Required Steps\n", "\n", "- Reading data from the catalog and appying the selection cuts.\n", "\n", "- Computing the tangential shear $e_{\\rm t}$ and inverse critical density $\\Sigma^{-1}_{\\rm crit}$. \n", "\n", "- $\\Delta \\Sigma (R)$ measurements using cKDTree and writing the output to a file.\n", "\n", "- Plotting the weak lensing signal.\n" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "#loading the required packages\n", "%matplotlib inline\n", "import pandas as pd\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.spatial import cKDTree\n", "from astropy.cosmology import FlatLambdaCDM\n", "import glob\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Reading data from the catalog and appying the selection cuts.\n", "\n", "We are going to use the selection cuts on lensing sample in the day-1 session. The below code can be use for other cuts too but for now lets stick to defaults as given below." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# selection cut on the lens sample\n", "def lens_select(zmin=0.1, zmax=0.33, lammin=55, lammax=100):\n", " #please check the file path properly \n", " data = pd.read_csv('/home/idies/workspace/Storage/divyar/IAGRG_2022/DataStore/redmapper.dat', delim_whitespace=1)\n", " #sample selection cut\n", " idx = (data['lambda']>lammin) & (data['lambda']<=lammax)\n", " idx = idx & (data['zred']>zmin) & (data['zred']<=zmax)\n", " ra = data['ra'].values[idx]\n", " dec = data['dec'].values[idx]\n", " zred = data['zred'].values[idx]\n", " #as we have no weights to apply we set them to unity\n", " wgt = ra*1.0/ra\n", " print('number of lenses=%d'%len(ra))\n", " return ra, dec, zred, wgt\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Similar we define a function for sources and collect the data from it.\n", "Please note that there are many columns in source file but for now we are only using some of them. \n", "Here we are only using \n", "\n", "- ra_gal, decgal : ra and dec for the sources\n", "- e1gal, e2gal: decribes the shapes of the sources\n", "- wgal, rms_egal: weights and Intrinsic shape dispersion per component\n", "- zphotgal: redshift of the sources\n", "\n", "For now we will neglect the data in other columns. As it is used for correcting the biases for our measurements and we will decribe how to use them at the end of this session. " ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "def read_sources(ifil):\n", " # various columns in sources \n", " # ragal, decgal, e1gal, e2gal, wgal, rms_egal, mgal, c1gal, c2gal, R2gal, zphotgal\n", " data = pd.read_csv(ifil, delim_whitespace=1).values\n", " zphotgal = data[:,-1]\n", " # sanity checks on the sources data\n", " idx = (np.sum(np.isnan(data), axis=1)==0) & (zphotgal>0)\n", " datagal = np.zeros((np.sum(idx),7))\n", " datagal[:,:6] = data[idx,:6]\n", " datagal[:,6] = data[idx,-1]\n", " # collects only - ragal, decgal, e1gal, e2gal, wgal, rms_egal, zphotgal\n", " return datagal" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Computing the tangential shear $e_{\\rm t}$ and inverse critical density $\\Sigma^{-1}_{\\rm crit}$ " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we will follow the formalism described in Prof Surhud More's lectures. We will first code up function to compute tangential component of the ellipticity given the positions for the lenses ($\\alpha_l$, $\\delta_l$) and sources($\\alpha_s$, $\\delta_s$) along with the shape measurements ($e_1$,$e_2$) for the sources.\n", "\n", "We will first compute angle $\\theta$ between a given lens-source pair.\n", "$$\\cos \\theta = \\cos \\delta_l \\cos \\delta_s \\cos(\\alpha_l - \\alpha_s) + \\sin\\delta_l \\sin \\delta_s $$\n", "\n", "where $\\delta_{l,s}$ and $\\alpha_{l,s}$ are ra and dec for lens(l) and source(s). The tangential component of ellipticity $e_t$ is given as\n", "\n", "$$e_t = - e_1 \\cos 2\\phi - e_2 \\sin 2\\phi$$\n", "\n", "$$ \\sin \\phi = \\frac{-\\sin \\delta_l \\cos \\delta_s + \\cos \\delta_l \\sin \\delta_s \\cos(\\alpha_s - \\alpha_l)}{|\\sin\\theta|} $$\n", "\n", "$$\\cos \\phi = \\frac{\\cos\\delta_l \\sin(\\alpha_s - \\alpha_l)}{|\\sin\\theta|}$$\n", "\n", "\n", "For more info related to above equations refer to Lecture Notes for [Tangential-shear-computation](https://surhudm.github.io/Weaklensing_IAGRG/Weak_gravitational_lensing.html#Tangential-shear-computation)\n", "\n", "\n", "Coding up these equations requires trigonometric functions from numpy package. Please check them out - [sin]( https://numpy.org/doc/stable/reference/generated/numpy.sin.html), [cos](https://numpy.org/doc/stable/reference/generated/numpy.cos.html#numpy.cos)\n", "\n", "**Please note that these functions by default uses angles in radians and here we are working with catalog data with (ra,dec) in degrees. So, we need to do the needful conversion first** " ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "# lra, ldec - lenses position\n", "# sra, sdec - sources position\n", "# se1 and se2 - source ellipticities\n", "def get_et(lra, ldec, sra, sdec, se1, se2):\n", " # angle_in_radian = angle_in_degrees * np.pi/180\n", " lra = ??\n", " # fill up the blanks\n", " ldec = ??\n", " sra = ??\n", " sdec = ?? \n", " # c_theta = cos_theta, s_theta = sin_theta\n", " # use the equations above and complete the expressions\n", " c_theta = np.cos(??)*np.cos(sdec) * ?? + np.sin(ldec)*np.sin(??)\n", " s_theta = np.sqrt(1-c_theta**2)\n", "\n", " # phi to get the compute the tangential shear\n", " c_phi = ?? *1.0/s_theta\n", " s_phi = (-np.sin(ldec)*np.cos(sdec) + ?? )*1.0/s_theta\n", " # tangential shear\n", " e_t = - se1*(??) - se2*(??)\n", "\n", " return e_t\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "Exercise: \n", " \n", "- Tell what are you getting the output for $e_t$ if you use input lra=0.0, ldec=0.0, sra=0.123, sdec=0.045, se1 = 4.5e-2, se2 = 1.7e-2.\n", "- running command : **print(get_et(lra=0.0, ldec=0.0, sra=0.123, sdec=0.045, se1 = 4.5e-2, se2 = 1.7e-2))** in the next cell below the function defination cell. Please add cell below using the insert option on the top. \n", " \n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "we will now move on to writing a function to get $\\Sigma^{-1}_{\\rm crit} (z_l,z_s)$ given the lense redshift ($z_l$) and source redshift ($z_s$) and it also needs a instance of astropy class (cc). Remember we need to create a instance for the astropy cosmo class given the input cosmological parameters. In the current code structure we will initiate the astopy cosmo class in our main code.\n", "\n", "To review you can refer to Day-1 docs page to see how astropy is used to get cosmological distances. \n", "Please note that here we are working in comoving coordinates to do the signal computations so the corresponding $\\Sigma^{-1}_{\\rm crit} (z_l,z_s)$ is given as\n", "\n", "$$\\Sigma^{-1}_{\\rm crit}(z_l,z_s) = \\frac{4\\pi G}{c^2} \\frac{d_{\\rm ang}(z_l) d_{\\rm ang}(z_l, z_s) (1 + z_l)^2} {d_{\\rm ang}(z_s)}$$\n", "\n", "where $d_{\\rm ang}(z)$ is the angular diameter distance for redshift $z$. \n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def get_sigma_crit_inv(lzred, szred, cc):\n", " # some important constants for the sigma crit computations\n", " gee = 4.301e-9 #km^2 Mpc M_sun^-1 s^-2 gravitational constant\n", " cee = 3e5 #km s^-1\n", " # sigma_crit_calculations for a given lense-source pair\n", " sigm_crit_inv = ?? * cc.angular_diameter_distance_z1z2(lzred, szred).value * (1.0 + lzred)**?? * 1.0/cc.angular_diameter_distance(szred).value\n", " sigm_crit_inv = sigm_crit_inv * 4*np.pi*gee*1.0/cee**2 \n", " sigm_crit_inv = 1e12*sigm_crit_inv #esd's are in pc not in Mpc\n", "\n", " return sigm_crit_inv\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "Exercise: \n", " \n", "- In the cell below the function defination, first initiate the cosmo class from astropy using \"**cc = FlatLambdaCDM(H0=100, Om0=0.999)**\".\n", "- running command : **print(get_sigma_crit_inv(lzred=0.33, szred=0.8, cc=cc))** \n", " \n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Remember how to convert from ra,dec (degrees) to x,y,z coordinates." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def get_xyz(ra, dec):\n", " ra = ra*np.pi/??\n", " dec = dec*??\n", " x = np.cos(dec)*np.cos(ra)\n", " y = np.cos(??)*np.sin(??)\n", " z = np.sin(dec) \n", " return x, y, z" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "Exercise: \n", " \n", "- Tell the output of : **print(get_xyz(30, 60))** \n", " \n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### $\\Delta \\Sigma$ measurements using cKDTree and writing the output to a file." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our aim in this section is to write the main function to measure weak lensing signal $\\Delta \\Sigma (R)$. As discussed in the lectures the $\\Delta \\Sigma (R)$ for a comoving projected radial bin $R$ as given by\n", "\n", "$$\\Delta \\Sigma (R) = \\frac{\\sum_{ls} w_{\\rm ls} e_{t,ls} (\\Sigma^{-1}_{\\rm crit})^{-1}}{2\\mathcal{R} \\sum_{ls} w_{\\rm ls}}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $w_{ls}$ is the source-lense weights given as $w_{ls} = w_l w_s \\Sigma^{-2}_{\\rm crit}$ and the $\\mathcal{R}$ is the responsivity given as $\\mathcal{R} = 1 - \\frac{\\sum_{\\rm ls} w_{ls} e^2_{\\rm rms}}{\\sum_{\\rm ls} w_{ls}}$. The $e_{\\rm rms}$ is the intrinsic shape dispersion given by the shape catalog for individual galaxies.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Decription of algorithm used in the code below (run_pipe) to get the signal $\\Delta \\Sigma$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1. Global Input Parameters:\n", " 1. **Omegam** : Omega matter - will be used in cosmological distance calculation, default = 0.315\n", " 2. **rmin** : Minimum projected radius - will be used for setting projected radial range, default = 0.2\n", " 3. **rmax** : Maximum projected radius - will be used for setting projected radial range, default = 2.0\n", " 4. **nbins** : number of projected radial bins - use to setup the logarithmic binning, default = 10\n", " 5. **zdiff** : a small offset added to lense redshift to select cleaner background sources $z_s > (z_l + {\\rm zdiff})$, default = 0.4 \n", " 6. **outputfile**: dat file name for the output, default = 'iagrg_dsigma.dat'\n", "\n", "2. **Initialize the astropy cosmo then getting projected radial binning array ready**\n", "\n", "3. We are computing the numerators and denominators for the signals separately and they have to be initialize first:\n", " 1. **sumdsig_num** : array for $\\Delta\\Sigma$ numerator\n", " 2. **sumdsigsq_num**: array for estimating the shape noise (covered in lectures) numerators\n", " 3. **sumwls** : array for summed lense-source weights\n", " 4. **sumwls_resp** : array for $\\Delta\\Sigma$ numerator(we can push the 1 into summed w_ls - look at the expression)\n", "\n", "4. Collect lenses data and convert the position (ra,dec) --> (x,y,z) using the above defined functions. Then use these x,y,z to get a lens tree using **cKDtree**.\n", "\n", "5. We also need to set a **maximum search radius** for the **query_ball_point** (refer to day-2 tutorials). The maximum is theratio between maximum projected radius rmax and minimum comoving distance to our lense sample - **rmax/(comoving_distance(zlmin))**, zlmin - minimum lense redshift. \n", "\n", "6. Now we use glob package to collect the file paths to our many many files. We loop over these path and read each file using the **read_source** function define above. \n", "\n", "7. Then we take the source ra, dec and convert them to x,y,z so that we can **query** around each source using lense tree and collect lenses within **maximum search radius**. This will give us a list of ids of lenses for each sources. \n", "\n", "8. Once we have a list of indices of lenses for all galaxies in the file. We then loop over the galaxies in individual file to compute the array of **comoving projected radial distance** for the each galaxy with it's lenses that we got in the last step. Here we also put a cut that for a given lense-source pair with the source redshift $z_s$ and lense redshift $z_l$ , $z_s > z_l +$ zdiff and we also skip sources which doesn't have any lense around them.\n", "\n", "9. The projected distance array from last step then used to look for corresponding bin number in our binning scheme as described in the Global parameters. Then we use the above equations along with the above defined functions to get the values of array defined in step 3.\n", "\n", "10. As said in step 6, we run step 7-9 for all the source files in our directory and then write the resultant output in a file named by the Global Parameters." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def run_pipe(Omegam=0.315, rmin=0.2, rmax=2.0, nbins=10, zdiff=0.4, outputfile = 'iagrg_dsigma.dat'):\n", " #set the cosmology with omegaM parameter \n", " cc = FlatLambdaCDM(H0=100, Om0=Omegam) # fixing H0=100 to set units in Mpc h-1\n", " \n", " # set the projected radial binning \n", " rmin = rmin\n", " rmax = rmax\n", " nbins = nbins #10 radial bins for our case\n", " rbins = np.logspace(np.log10(rmin), np.log10(rmax), nbins + 1)\n", " rdiff = np.log10(rbins[1]*1.0/rbins[0]) # difference in log binning\n", " \n", " # initializing arrays for signal compuations\n", " sumdsig_num = np.zeros(len(rbins[:-1]))\n", " sumdsigsq_num = np.zeros(len(rbins[:-1]))\n", " sumwls = np.zeros(len(rbins[:-1]))\n", " sumwls_resp = np.zeros(len(rbins[:-1]))\n", "\n", " # getting the lenses data\n", " lra, ldec, lred, lwgt = ?? \n", " # convert lense ra and dec into x,y,z cartesian coordinates\n", " lx, ly, lz = ?? \n", " # putting kd tree around the lenses\n", " lens_tree = ?? \n", " print('lenses tree is ready\\n')\n", " \n", " # setting maximum search radius\n", " dcommin = cc.comoving_distance(np.min(lred)).value\n", " dismax = (rmax*1.0/(dcommin)) \n", "\n", " # lets first catch the file list for sources\n", " sflist = np.sort(glob.glob('/home/idies/workspace/Storage/divyar/IAGRG_2022/DataStore/hsc/*.dat'))\n", "\n", " # Ready to pounce on the source data\n", " for ifil in sflist:\n", " # source data in datagal array using read source function\n", " datagal = ??\n", " Ngal = len(datagal[:,0]) # total number of galaxies in the source file\n", " \n", " # first two entries are ra and dec for the sources\n", " allragal = datagal[:,0]\n", " alldecgal = datagal[:,1]\n", " \n", " # ra and dec to x,y,z for sources\n", " allsx, allsy, allsz = ??\n", " \n", " # query in a ball around individual sources and collect the lenses ids with a maximum radius\n", " slidx = lens_tree.query_ball_point(np.transpose([allsx, allsy, allsz]), dismax) \n", " \n", " # looping over all the galaxies\n", " for igal in range(Ngal): \n", " ragal = datagal[igal,0]\n", " decgal = datagal[igal,1]\n", " e1gal = datagal[igal,2]\n", " e2gal = datagal[igal,3]\n", " wgal = datagal[igal,4]\n", " rms_egal = datagal[igal,5]\n", " zphotgal = datagal[igal,6]\n", " \n", " # array of lenses indices\n", " lidx = np.array(slidx[igal])\n", " \n", " # skip sources which doesn't have any lenses around them \n", " if len(lidx)==0:\n", " continue\n", " \n", " # selecting a cleaner background\n", " zcut = (lred[lidx] < (zphotgal - zdiff)) #only taking the foreground lenses\n", " \n", " # again skipping the onces which doesn't satisfy the above criteria\n", " if np.sum(zcut)==0.0:\n", " continue\n", " \n", " # collecting the data of lenses around individual source\n", " lidx = lidx[zcut] # this will catch the array indices for our lenses\n", " sra = ragal\n", " sdec = decgal\n", " \n", " l_ra = lra[lidx]\n", " l_dec = ldec[lidx]\n", " l_zred = lred[lidx] \n", " l_wgt = lwgt[lidx] \n", " \n", " sx, sy, sz = ?? # individual galaxy sra,sdec-->sx,sy,sz\n", " lx, ly, lz = ?? # individual galaxy lra,ldec-->lx,ly,lz\n", " \n", " # getting the radial separations for a lense source pair \n", " sl_sep = np.sqrt((lx - sx)**2 + (ly - sy)**2 + (lz - sz)**2)\n", " sl_sep = sl_sep * cc.comoving_distance(l_zred).value\n", " \n", " # here we are binning the separation in projected radial bins\n", " # ll will enumerate the lenses\n", " for ll,sep in enumerate(sl_sep):\n", " \n", " #check whether separation sits inside our projected radial range\n", " if seprmax: \n", " continue\n", " \n", " # finding the corresponding radial bins for separations \n", " rb = int(np.log10(sep*1.0/rmin)*1/rdiff)\n", " \n", " # get tangantial components given positions and shapes\n", " e_t = ??(lra = l_ra[ll], ldec = l_dec[ll], sra = ??, sdec = ??, se1 = ??, se2 = ??)\n", "\n", " # sigma_crit_calculations for a given lense-source pair with cc as astropy cosmo instance\n", " sigm_crit_inv = ??(l_zred[ll], ??, cc)\n", "\n", " # following equations given in the surhud's lectures \n", " w_ls = l_wgt[ll] * wgal * (??)**??\n", " w_ls_by_av_sigc_inv = l_wgt[ll] * ?? * sigm_crit_inv\n", "\n", " # separate numerator and denominator computation \n", " sumdsig_num[rb] += w_ls_by_av_sigc_inv * e_t\n", " sumdsigsq_num[rb] += (w_ls_by_av_sigc_inv * e_t)**2\n", " sumwls[rb] += w_ls\n", " sumwls_resp[rb] += w_ls * (1-rms_egal**2)\n", "\n", " print(ifil)\n", " \n", " outputfile = outputfile \n", " fout = open(outputfile, \"w\")\n", " fout.write(\"# 0:rmin/2+rmax/2 1:DeltaSigma 2:SN_ErrDeltaSigma\\n\")\n", " for i in range(len(rbins[:-1])):\n", " rrmin = rbins[i]\n", " rrmax = rbins[i+1]\n", " Resp = sumwls_resp[i]*1.0/sumwls[i]\n", " dsigma = sumdsig_num[i]*1.0/sumwls[i]/2./Resp # follow from the equation above\n", " dsigma_err = np.sqrt(sumdsigsq_num[i])*1.0/sumwls[i]/2./Resp\n", " \n", " fout.write(\"%le\\t%le\\t%le\\n\"%(rrmin/2.0+rrmax/2.0, dsigma, dsigma_err))\n", " fout.write(\"#OK\") \n", " fout.close()\n", " \n", " return 0" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "run_pipe()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plotting the weak lensing signal" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "dat = np.loadtxt('iagrg_dsigma.dat')\n", "\n", "plt.errorbar(dat[:,0], dat[:,1], yerr=dat[:,2], fmt='.', capsize=3, label='Data')\n", "plt.legend()\n", "\n", "plt.xlabel(r'$R[{\\rm h^{-1}Mpc}]$')\n", "plt.ylabel(r'$\\Delta\\Sigma [{\\rm h M_\\odot pc^{-2}}]$')\n", "plt.xscale('log')\n", "plt.yscale('log')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The bias corrected measurements " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For doing a bias corrected measurements one has to add few more things in our $\\Delta \\Sigma$ expression. Interested reader can check this expression for $\\Delta \\Sigma$ in [arXiv:2107.05641](https://arxiv.org/abs/2107.05641) and can modify the above code to get the bias corrected measurements.\n", "\n", "Please note that within each source file we have already provided you the data coloumns required to do these computations. \n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "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.7.3" } }, "nbformat": 4, "nbformat_minor": 4 }