CMB_School_Part_06.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"https://colab.research.google.com/github/jeffmcm1977/CMBAnalysis_SummerSchool/blob/master/CMB_School_Part_06.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Map Analysis II\n",
    "\n",
    "###  Renée Hložek and Jeff McMahon\n",
    "\n",
    "## Using our tools on real data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this self study exercise, we are now going to use some public data from the ACT collaboration, compute power spectra and compare them to theory power spectra."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!python -c \"import cmb_modules\" || ( \\\n",
    "    wget https://github.com/jeffmcm1977/CMBAnalysis_SummerSchool/raw/master/cmb_school.tar.gz && \\\n",
    "    tar xzvf cmb_school.tar.gz \\\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib\n",
    "import sys\n",
    "import matplotlib.cm as cm\n",
    "import matplotlib.mlab as mlab\n",
    "import matplotlib.pyplot as plt\n",
    "import astropy.io.fits as fits\n",
    "\n",
    "import constants as cs # the constants module\n",
    "import cmb_modules # the module of functions\n",
    "\n",
    "N = cs.N\n",
    "c_min = cs.c_min\n",
    "c_max = cs.c_max\n",
    "X_width =cs.X_width\n",
    "Y_width = cs.Y_width\n",
    "beam_size_fwhp = cs.beam_size_fwhp\n",
    "\n",
    "pix_size = cs.pix_size\n",
    "\n",
    "Number_of_Sources  = cs.Number_of_Sources\n",
    "Amplitude_of_Sources = cs.Amplitude_of_Sources\n",
    "Number_of_Sources_EX = cs.Number_of_Sources_EX\n",
    "Amplitude_of_Sources_EX = cs.Amplitude_of_Sources_EX\n",
    "\n",
    "Number_of_SZ_Clusters  = cs.Number_of_SZ_Clusters\n",
    "Mean_Amplitude_of_SZ_Clusters = cs.Mean_Amplitude_of_SZ_Clusters\n",
    "SZ_beta = cs.SZ_beta\n",
    "SZ_Theta_core = cs.SZ_Theta_core\n",
    "\n",
    "white_noise_level = cs.white_noise_level\n",
    "atmospheric_noise_level = cs.atmospheric_noise_level\n",
    "one_over_f_noise_level = cs.one_over_f_noise_level"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reading in data and making simulated maps\n",
    "Let's start by reading in a map and computing the power spectrum of the map, and then comparing to a theory map of a similar size/shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map = fits.getdata('ACT_Season3_148_GHz_submap.fits')\n",
    "header = fits.getheader('ACT_Season3_148_GHz_submap.fits')\n",
    "ra0 = 310.006000\n",
    "ra1 = 360.001000\n",
    "dec0 = -1.996904\n",
    "dec1 = 1.988649\n",
    "Nx= header['NAXIS1']\n",
    "Ny= header['NAXIS2']\n",
    "\n",
    "pix_x = 60.*(ra1-ra0)/Nx\n",
    "pix_y = 60.*(dec1-dec0)/Ny\n",
    "print(pix_x, pix_y) # x and y pixel scales in arcmins\n",
    "plt.figure(figsize=(10,10))\n",
    "p = cmb_modules.Plot_CMB_Map(map,c_min,c_max,Nx,Ny)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see that these ACT maps are long stripes, so we will cut a square patch out of the long strip and compute the PS on that patch."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N=int(Ny) # take only one dimension\n",
    "map_patch = map[0:int(N), 0:int(N)] #trimming the map down to size\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Making a window function, taking the shorter dimension for now\n",
    "window = (cosine_window(N))\n",
    "appodized_map = window * map_patch\n",
    "p=cmb_modules.Plot_CMB_Map(appodized_map,c_min,c_max,N, N)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### We can see that there is a lot of power on large scales from the atmosphere. What do we expect to see in the power spectrum?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Make a CMB map\n",
    "ell, DlTT = np.loadtxt(\"CAMB_fiducial_cosmo_scalCls.dat\", usecols=(0, 1), unpack=True) \n",
    "\n",
    "CMB_T = cmb_modules.make_CMB_T_map(N,pix_size,ell,DlTT)\n",
    "\n",
    "## make a point source map\n",
    "PSMap = cmb_modules.Poisson_source_component(N,pix_size,Number_of_Sources,Amplitude_of_Sources) \n",
    "PSMap += cmb_modules.Exponential_source_component(N,pix_size,Number_of_Sources_EX,Amplitude_of_Sources_EX)\n",
    "\n",
    "## make an SZ map\n",
    "SZMap,SZCat = cmb_modules.SZ_source_component(N,pix_size,Number_of_SZ_Clusters,Mean_Amplitude_of_SZ_Clusters,SZ_beta,SZ_Theta_core,False)\n",
    "\n",
    "## add them all together to get the sky map at a single freuqency\n",
    "total_map = CMB_T + PSMap + SZMap\n",
    "\n",
    "## incorperate the impact of the instrument\n",
    "    ## beam\n",
    "CMB_T_convolved = cmb_modules.convolve_map_with_gaussian_beam(N,pix_size,beam_size_fwhp,total_map)\n",
    "    ## noise\n",
    "Noise = cmb_modules.make_noise_map(N,pix_size,white_noise_level,atmospheric_noise_level,one_over_f_noise_level)\n",
    "\n",
    "total_map_plus_noise = CMB_T_convolved + Noise\n",
    "\n",
    "## plot the result\n",
    "#p = cmb_modules.Plot_CMB_Map(total_map_plus_noise,c_min,c_max,N,N)\n",
    "\n",
    "appodized_theory_map = window * total_map_plus_noise\n",
    "\n",
    "p = cmb_modules.Plot_CMB_Map(appodized_theory_map,c_min,c_max,N, N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## We will now take the power spectrum of this apodised map, and the theory version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#### parameters for setting up the spectrum\n",
    "delta_ell = 10# cs.delta_ell\n",
    "ell_max = cs.ell_max\n",
    "#ell, DlTT = np.loadtxt(\"CAMB_fiducial_cosmo_scalCls.dat\", usecols=(0, 1), unpack=True) \n",
    "#plt.clf()\n",
    "\n",
    "if max(ell)< ell_max: \n",
    "        print('WARNING: Your theory curves end before the binned ell_max')\n",
    "        \n",
    "binned_ell_theory, binned_spectrum_theory = \\\n",
    "    cmb_modules.calculate_2d_spectrum(appodized_theory_map,delta_ell,ell_max,pix_size,N)\n",
    "binned_ell_dat, binned_spectrum_dat = \\\n",
    "    cmb_modules.calculate_2d_spectrum(appodized_map,delta_ell,ell_max,pix_size,N)\n",
    "\n",
    "#print binned_ell, binned_spectrum\n",
    "#print np.mean(binned_ell), np.shape(binned_ell), np.max(binned_ell), np.min(binned_ell)\n",
    "plt.semilogy(binned_ell_theory,binned_spectrum_theory*\\\n",
    "                 binned_ell_theory * (binned_ell_theory+1.)/2. / np.pi, label='theory')\n",
    "plt.semilogy(binned_ell_dat,binned_spectrum_dat*\\\n",
    "                binned_ell_dat * (binned_ell_dat+1.)/2. / np.pi, label='data')\n",
    "#plt.clf()\n",
    "leg = plt.legend()\n",
    "leg.draw_frame(False)\n",
    "plt.semilogy(ell,DlTT)\n",
    "plt.ylabel('$D_{\\ell}$ [$\\mu$K$^2$]')\n",
    "plt.xlabel('$\\ell$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "source": [
    "Notice that the power spectrum is biased high on large scales - as we expected it to be.\n",
    "How might you remove that large scale power before computing the power spectrum?\n",
    "<font color='red'>EXERCISE: </font>  Write code to filter out the large scale modes ell < ell_min before computing the power spectrum and show the spectra for a few different values for ell_min. Also, can you think about how to remove the ringing the power on small scales? Discuss the shape of the apodising window and the large scale power leaking to small scales."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Your code and plots here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>EXERCISE: </font>  We cut out a particular patch from the ACT strip. Cut out different patches and view them. By using the power spectrum code in the modules file, compute both the auto and the cross spectra for this map. Discuss how you would estimate the error bars on this power spectrum, by comparing the auto with the cross power spectra."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>EXERCISE: </font>  Given the value of the PS shown here on small scales, we can estimate the noise level in the ACT maps for this season. Plot noise curves for a few values and show here, and do a simple \"chi-by-eye\" fit for the noise level."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Your discussion here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>EXERCISE: </font>  Plot the 2D power spectrum, both of your input theory+ noise and the ACT data. What does that tell you about real-world noise from CMB Experiments?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Frequency dependence of maps and foregrounds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are going to look at ACT maps on the same patch of the sky, but in a different frequency band."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map220 = fits.getdata('ACT_Season3_220_GHz_submap.fits')\n",
    "header = fits.getheader('ACT_Season3_220_GHz_submap.fits')\n",
    "ra0 = 310.006000\n",
    "ra1 = 360.001000\n",
    "dec0 = -1.996904\n",
    "dec1 = 1.988649\n",
    "Nx= header['NAXIS1']\n",
    "Ny= header['NAXIS2']\n",
    "\n",
    "pix_x = 60.*(ra1-ra0)/Nx\n",
    "pix_y = 60.*(dec1-dec0)/Ny\n",
    "print(pix_x, pix_y) # x and y pixel scales in arcmins\n",
    "plt.figure(figsize=(10,10))\n",
    "p = cmb_modules.Plot_CMB_Map(map220,c_min,c_max,Nx,Ny)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N=int(Ny) # take only one dimension\n",
    "map_patch220 = map220[0:int(N), 0:int(N)] #trimming the map down to size\n",
    "window = (cosine_window(N))\n",
    "appodized_map220 = window * map_patch220\n",
    "\n",
    "binned_ell_dat220, binned_spectrum_dat220 = cmb_modules.calculate_2d_spectrum(appodized_map220,delta_ell,ell_max,pix_size,N)\n",
    "plt.semilogy(binned_ell_dat,binned_spectrum_dat* binned_ell_dat * (binned_ell_dat+1.)/2. / np.pi, label='148')\n",
    "plt.semilogy(binned_ell_dat220,binned_spectrum_dat220* binned_ell_dat220 * (binned_ell_dat220+1.)/2. / np.pi, label='220')\n",
    "#plt.clf()\n",
    "leg = plt.legend()\n",
    "leg.draw_frame(False)\n",
    "plt.semilogy(ell,DlTT)\n",
    "plt.ylabel('$D_{\\ell}$ [$\\mu$K$^2$]')\n",
    "plt.xlabel('$\\ell$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>EXERCISE: </font> The spectrum is biased high from point sources on all scales. Use the code from the previous notebook to determine the noise bias for this spectrum and remove it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "source": [
    "<font color='red'>EXERCISE: </font> Apply the techniques from Part Five and search for point sources and SZ clusters in the map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}