stack.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Author: M. Riley Owens (GitHub: mrileyowens)\n",
    "\n",
    "# This file makes stacks of the spectra of the LyC-leaking \n",
    "# region and the non-LyC-leaking regions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import glob\n",
    "\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def stack():\n",
    "    \n",
    "    '''\n",
    "    Stacks individual MagE spectra into two stacked spectra: one of the LyC-leaking region, and another \n",
    "    of the non-LyC-leaking regions\n",
    "    '''\n",
    "\n",
    "    def get_data(file):\n",
    "\n",
    "        '''\n",
    "        Extracts the wavelength, flux density, and uncertainty from a given MagE spectrum\n",
    "\n",
    "        Parameters:\n",
    "            file : str\n",
    "                Complete file path to the MagE spectrum\n",
    "\n",
    "        Returns:\n",
    "            w : numpy.ndarray\n",
    "                Observed wavelength bins of the spectrum (angstroms)\n",
    "            f : numpy.ndarray\n",
    "                Observed flux densities of the spectrum (erg/s/cm^2/angstrom)\n",
    "            n : numpy.ndarray\n",
    "                Gaussian 1σ observed flux density uncertainties of the spectrum\n",
    "        '''\n",
    "\n",
    "        # Get the observed wavelengths (angstroms), flux densities (erg/s/cm^2/Hz), and Gaussian 1σ flux density uncertainties of the spectrum\n",
    "        w, f, n = np.loadtxt(file, delimiter='\\t', comments=('#', 'w'), usecols=(0,1,2), unpack=True)\n",
    "\n",
    "        # Only keep array entries where the flux densities and uncertainties are physically sensible\n",
    "        w = w[f < 1e-20]\n",
    "        n = n[f < 1e-20]\n",
    "        f = f[f < 1e-20]\n",
    "\n",
    "        w = w[n < 1e-20]\n",
    "        f = f[n < 1e-20]\n",
    "        n = n[n < 1e-20]\n",
    "\n",
    "        # Convert the flux densities and uncertainties to wavelength space\n",
    "        f = f * 2.998e18 / np.square(w)\n",
    "        n = n * 2.998e18 / np.square(w)\n",
    "\n",
    "        return w, f, n\n",
    "\n",
    "    # Establish basic directories\n",
    "    home = os.getcwd()\n",
    "    data = f'{home}/data'\n",
    "\n",
    "    # Get the file paths of the MagE spectra. This line is not optimal, since it will include\n",
    "    # the stacked spectra if the notebook executes again after creating them, which could\n",
    "    # cause problems later in the notebook. I (M. Riley Owens) am leaving it as is with this \n",
    "    # danger noted\n",
    "    files = glob.glob(f'{data}/mage/*.txt')\n",
    "\n",
    "    # Dictionary of the spectroscopic redshifts of each spectrum, as \n",
    "    # measured by Mainali et al. (2022) (ApJ, 940, 160) from the narrow \n",
    "    # component of the [O III] 5007 Angstrom nebular emission line.\n",
    "    z_dict = {\n",
    "        'sunburst_M-0-comb1_MWdr' : 2.37014,\n",
    "        'sunburst_M-2-comb1_MWdr' : 2.37017,\n",
    "        'sunburst_M-3-comb1_MWdr' : 2.37025,\n",
    "        'sunburst_M-4-comb1_MWdr' : 2.37073,\n",
    "        'sunburst_M-5-comb1_MWdr' : 2.37086,\n",
    "        'sunburst_M-6-comb1_MWdr' : 2.37021,\n",
    "        'sunburst_M-7-comb1_MWdr' : 2.37044,\n",
    "        'sunburst_M-8-comb1_MWdr' : 2.37024,\n",
    "        'sunburst_M-9-comb1_MWdr' : 2.37030\n",
    "    }\n",
    "\n",
    "    # Make empty arrays for the final stacked spectra,\n",
    "    # which will be filled in later\n",
    "    w_leaker = np.array([])\n",
    "    f_leaker = np.array([])\n",
    "    n_leaker = np.array([])\n",
    "\n",
    "    w_nonleaker = np.array([])\n",
    "    f_nonleaker = np.array([])\n",
    "    n_nonleaker = np.array([])\n",
    "\n",
    "    # For each individual spectrum\n",
    "    for i, file in enumerate(files):\n",
    "\n",
    "        # If the spectrum is one of the LyC leaker (excluding slit M0\n",
    "        # because of poor observing conditions leading to a poor fluxing)\n",
    "        if any(string in file for string in ['M-2','M-7','M-8','M-9']):\n",
    "\n",
    "            # Get the observed wavelength (angstroms), flux densities (erg/s/cm^2/angstrom), \n",
    "            # and Gaussian 1σ flux density uncertainties of the MagE spectrum\n",
    "            w, f, n = get_data(file)      \n",
    "\n",
    "            # Place the wavelength bins in the rest frame\n",
    "            w = w / (1 + z_dict[os.path.basename(file)[:-4]])\n",
    "\n",
    "            # Normalize the flux densities and uncertainties by the median flux density\n",
    "            # between rest-frame 1267 - 1276 angstroms, just like the continuum fits to\n",
    "            # the individual MagE spectra created by John Chisholm\n",
    "            n = n / np.median(f[(w >= 1267) & (w <= 1276)])\n",
    "            f = f / np.median(f[(w >= 1267) & (w <= 1276)])\n",
    "\n",
    "            # If the final arrays of the stacked LyC leaker spectrum have not\n",
    "            # been instantiated yet, do so now\n",
    "            if len(w_leaker) == 0:\n",
    "\n",
    "                w_leaker = w\n",
    "                f_leaker = f\n",
    "                n_leaker = n\n",
    "\n",
    "            else:\n",
    "\n",
    "                # Interpolate the data at the wavelength bins of the stacked spectrum\n",
    "                f = np.interp(w_leaker, w, f, left=0, right=0)\n",
    "                n = np.interp(w_leaker, w, n, left=0, right=0)\n",
    "\n",
    "                # Add the flux density and flux density uncertainties to the stacked spectrum\n",
    "                # appropriately\n",
    "                f_leaker += f\n",
    "                n_leaker = np.sqrt((n_leaker)**2 + (n)**2)\n",
    "\n",
    "        # If the spectrum is of a non-LyC-leaking region (excluding slit M3 because it captures a highly compact object; see Vanzella \n",
    "        # et al. (2020) (MNRAS, 499, L67), Diego et al. (2022) (A&A, 665, A134), or Sharon et al. (2022) (ApJ, 941, 203))\n",
    "        elif any(string in file for string in ['M-5','M-4','M-6']):\n",
    "\n",
    "            # Get the observed wavelength (angstroms), flux densities (erg/s/cm^2/angstrom), \n",
    "            # and Gaussian 1σ flux density uncertainties of the MagE spectrum\n",
    "            w, f, n = get_data(file)\n",
    "\n",
    "            # Place the wavelength bins in the rest frame\n",
    "            w = w / (1 + z_dict[os.path.basename(file)[:-4]])\n",
    "\n",
    "            # Normalize the flux densities and uncertainties by the median flux density\n",
    "            # between rest-frame 1267 - 1276 angstroms, just like the continuum fits to\n",
    "            # the individual MagE spectra created by John Chisholm\n",
    "            n = n / np.median(f[(w >= 1267) & (w <= 1276)])\n",
    "            f = f / np.median(f[(w >= 1267) & (w <= 1276)])\n",
    "\n",
    "            # If the final arrays of the stacked non-LyC-leaking spectrum have not\n",
    "            # been instantiated yet, do so now\n",
    "            if len(w_nonleaker) == 0:\n",
    "\n",
    "                w_nonleaker = w\n",
    "                f_nonleaker = f\n",
    "                n_nonleaker = n\n",
    "\n",
    "            else:\n",
    "\n",
    "                # Interpolate the data at the wavelength bins of the stacked spectrum\n",
    "                f = np.interp(w_nonleaker, w, f, left=0, right=0)\n",
    "                n = np.interp(w_nonleaker, w, n, left=0, right=0)\n",
    "\n",
    "                # Add the flux density and flux density uncertainties to the stacked spectrum\n",
    "                # appropriately\n",
    "                f_nonleaker += f\n",
    "                n_nonleaker = np.sqrt((n_nonleaker)**2 + (n)**2)\n",
    "\n",
    "    # Normalize the two stacks by the number of composite spectra\n",
    "    f_leaker = f_leaker / 4\n",
    "    n_leaker = n_leaker / 4\n",
    "\n",
    "    f_nonleaker = f_nonleaker / 3\n",
    "    n_nonleaker = n_nonleaker / 3\n",
    "\n",
    "    # Save the stacked spectra to .txt files\n",
    "    np.savetxt(f'{data}/mage/sunburst_arc_leaker_stack_mage.txt', np.array([w_leaker, f_leaker, n_leaker]).T, delimiter='\\t', encoding='utf-8',\n",
    "        header=f'Stacked spectrum of the LyC-leaking region of the Sunburst Arc from the MagE slits M2, M7, M8, M9, created by stack.ipynb.' + '\\n' \n",
    "        + 'The stack excludes slit M0 because of poor observing conditions that caused a poor fluxing.' + '\\n' + '\\n'\n",
    "        + 'The columns below are, from left to right: rest wavelength (angstroms), flux density (arbitrary units in wavelength space), and the 1σ Gaussian uncertainty of the flux density.')\n",
    "    np.savetxt(f'{data}/mage/sunburst_arc_nonleaker_stack_mage.txt', np.array([w_nonleaker, f_nonleaker, n_nonleaker]).T, delimiter='\\t', encoding='utf-8',\n",
    "        header=f'Stacked spectrum of the non-LyC-leaking regions of the Sunburst Arc from the MagE slits M5, M4, M6, created by stack.ipynb.' + '\\n' \n",
    "        + 'The stack excludes slit M3 because it captures a highly compact region of the galaxy which may not be representative of it\\'s surrounding environment.' + '\\n'\n",
    "        + 'Consult Vanzella et al. (2020) (MNRAS:Letters, 499, 1, L67), Diego et al. (2022) (A&A, 665, A134), or Sharon et al. (2022) (ApJ, 941, 2, 203) for more information about the compact region.' + '\\n' + '\\n'\n",
    "        + 'The columns below are, from left to right: rest wavelength (angstroms), flux density (arbitrary units in wavelength space), and the 1σ Gaussian uncertainty of the flux density.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "stack()"
   ]
  }
 ],
 "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.12.2"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}