imaging_mode/imaging_mode_stage_1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='top'></a>\n",
    "# Imaging Mode Data Calibration: Part 1 - Ramps to Slopes\n",
    "---\n",
    "**Author**: Bryan Hilbert (hilbert@stsci.edu)| **Latest Update**: 7 May 2021"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    <h3><u><b>Notebook Goals</b></u></h3>\n",
    "    <ul>Working with the Stage 1 Calibration Pipeline, we will:</ul>    \n",
    "<ul>\n",
    "    <li>Look at the different ways to call the pipeline</li>\n",
    "    <li>Examine exactly what each pipeline step does to the science data</li>    \n",
    "</ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table of Contents\n",
    "* [Introduction](#intro)\n",
    "* [Pipeline Resources and Documentation](#resources)\n",
    "   * [Installation](#installation)\n",
    "   * [Reference Files](#reference_files)\n",
    "* [Imports](#imports)\n",
    "* [Convenience Functions](#convenience_functions)\n",
    "* [Download Data](#download_data)\n",
    "* [Methods for calling steps/pipelines](#calling_methods)\n",
    "* [Parameter reference files](#parameter_reffiles)\n",
    "* [calwebb_detector1 - Ramps to slopes](#detector1) \n",
    "   * [Run the entire pipeline](#detector1_at_once)\n",
    "       * [run() method](#run_method)\n",
    "       * [call() method](#call_method)\n",
    "       * [command line](#command_line)\n",
    "       * [Exercise 1](#exercise1)\n",
    "   * [Run the individual pipeline steps](#detector1_step_by_step)\n",
    "       * [The `Data Quality Initialization` step](#dq_init)\n",
    "       * [The `Saturation Flagging` step](#saturation)\n",
    "       * [The `Superbias Subtraction` step](#superbias)\n",
    "       * [The `Reference Pixel Subtraction` step](#refpix)\n",
    "           * [Exercise 2](#exercise2)\n",
    "       * [The `Linearity Correction` step](#linearity)\n",
    "       * [The `Persistence Correction` step](#persistence)\n",
    "       * [The `Dark Current Subtraction` step](#dc)\n",
    "       * [The `Cosmic Ray Flagging` step](#jump)\n",
    "       * [The `Ramp_Fitting` step](#ramp_fitting)\n",
    "* [Bonus Topic: Logging](#logging)\n",
    "* [Exercise solutions](#exercise_solutions)\n",
    "    * [Exercise 1 solution](#exercise1_solution)\n",
    "    * [Exercise 2 solution](#exercise2_solution)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='intro'></a>\n",
    "## Introduction\n",
    "\n",
    "This notebook covers part 1 of the imaging mode data calibration module. In this notebook we'll review Stage 1 of the JWST calibration pipeline, also known as *calwebb_detector1*. \n",
    "\n",
    "The [Stage 1 pipeline](https://jwst-pipeline.readthedocs.io/en/stable/jwst/pipeline/calwebb_detector1.html#calwebb-detector1) applies basic detector-level corrections to all exposure types (imaging, spectroscopic, coronagraphic, etc.). It is applied to one exposure at a time, beginning with an uncalibrated multiaccum ramp (*_uncal.fits file*). It is sometimes referred to as “ramps-to-slopes” processing. \n",
    "\n",
    "Each input raw data file is composed of one or more ramps (integrations) containing increasing count values from the non-destructive detector readouts. For details on multiaccum files and data collection, see the JDox page on [how up-the-ramp readouts work](https://jwst-docs.stsci.edu/understanding-exposure-times#UnderstandingExposureTimes-uptherampHowup-the-rampreadoutswork).\n",
    "\n",
    "The output is a corrected but uncalibrated countrate or slope image (*_rate.fits and _rateints.fits file*). In this case, \"calibrated/uncalibrated\" is based on the units of the data. Any data in units of DN or DN/sec are considered uncalibrated. After the flux calibration step of the Stage 2 pipeline is run and the data units become physical units (e.g. MJy/str), then the data are said to be calibrated.\n",
    "\n",
    "To illustrate how the steps of the pipeline change the input data, we will download a sample file and run it through the pipeline, examining the results at several places along the way.\n",
    "\n",
    "All JWST data, regardless of instrument and observing mode, are processed through the Stage 1 pipeline. The corrections performed are the same across all near-IR instruments. There are several additional MIRI-specific steps. For the purposes of this notebook, our example file will be NIRCam data. We will also provide an example MIRI file that can be used in a separate exercise."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='resources'></a>\n",
    "## Pipeline Resources and Documentation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are several different places to find information on installing and running the pipeline. This notebook will summarize each Stage 1 pipeline step. To find more in-depth instructions use the links below.\n",
    "\n",
    "* [JWST Documentation (JDox) for the Stage 1 pipeline](https://jwst-docs.stsci.edu/jwst-data-reduction-pipeline/algorithm-documentation/stages-of-processing/calwebb_detector1) including short a short summary of what each step does.\n",
    "\n",
    "* [High-level description of all pipeline stages and steps](https://jwst-pipeline.readthedocs.io/en/latest/jwst/pipeline/main.html) from the pipeline software documentation website.\n",
    "\n",
    "* [`jwst` package documentation](https://jwst-pipeline.readthedocs.io/en/latest/jwst/introduction.html) including how to run the pipeline, input/output files, etc.\n",
    "\n",
    "* [`jwst` package GitHub repository, with installation instructions](https://github.com/spacetelescope/jwst/blob/master/README.md)\n",
    "\n",
    "* [**Help Desk**](https://stsci.service-now.com/jwst?id=sc_cat_item&sys_id=27a8af2fdbf2220033b55dd5ce9619cd&sysparm_category=e15706fc0a0a0aa7007fc21e1ab70c2f): **If you have any questions or problems regarding the pipeline, submit a ticket to the Help Desk**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='installation'></a>\n",
    "### Installation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "    During the JWebbinar, we will be working in a pre-existing environment where the <b>jwst</b> package has already been installed, so you won't need to install it yourself.\n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-block alert-warning\">\n",
    "If you wish to run this notebook outside of this JWebbinar, you will have to first install the <b>jwst</b> package.<br>\n",
    "\n",
    "For more detailed instructions on the various ways to install the package, see the [installation instructions](https://github.com/spacetelescope/jwst/blob/master/README.md) on GitHub.\n",
    "\n",
    "The easiest way to install the pipeline is via `pip`. Below we show how to create a new conda environment, activate that environment, and then install the latest released version of the pipeline. You can name your environment anything you like. In the lines below, replace `<env_name>` with your chosen environment name.\n",
    "\n",
    ">`conda create -n <env_name> python`<br>\n",
    ">`conda activate <env_name>`<br>\n",
    ">`pip install jwst`\n",
    "\n",
    "If you wish to install the development version of the pipeline, which is more recent than (but not as well tested compared to) the latest released version:\n",
    "\n",
    ">`conda create -n <env_name> python`<br>\n",
    ">`conda activate <env_name>`<br>\n",
    ">`pip install git+https://github.com/spacetelescope/jwst`\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='reference_files'></a>\n",
    "### Reference Files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Calibration reference files](https://jwst-docs.stsci.edu/data-processing-and-calibration-files/calibration-reference-files) are a collection of FITS and ASDF files that are used to remove instrumental signatures and calibrate JWST data. For example, the dark current reference file contains a multiaccum ramp of dark current signal to be subtracted from the data during the dark current subtraction step. \n",
    "\n",
    "When running a pipeline or pipeline step, the pipeline will automatically look for any required reference files in a pre-defined local directory. If the required reference files are not present, they will automatically be downloaded from the Calibration Reference Data System (CRDS) at STScI.\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "    During the JWebbinar, our pre-existing existing environment is set up to correctly use and store calibration reference files, and you do not need to set the environment variables below.\n",
    "</div>\n",
    "    \n",
    "<div class=\"alert alert-block alert-warning\">\n",
    "If you wish to run this notebook outside of this JWebbinar, you will have to specify a local directory in which to store reference files, along with the server to use to download the reference files from CRDS. To accomplish this, there are two environment variables that should be set prior to calling the pipeline. These are the CRDS_PATH and CRDS_SERVER_URL variables. In the example below, reference files will be downloaded to the \"crds_cache\" directory under the home directory.\n",
    "\n",
    ">`$ export CRDS_PATH=$HOME/crds_cache`<br>\n",
    ">`$ export CRDS_SERVER_URL=https://jwst-crds.stsci.edu`<br>\n",
    "OR:<br>\n",
    "`os.environ[\"CRDS_PATH\"] = \"/user/myself/crds_cache\"`<br>\n",
    "`os.environ[\"CRDS_SERVER_URL\"] = \"https://jwst-crds.stsci.edu\"`<br>\n",
    "\n",
    "The first time you run the pipeline, the CRDS server should download all of the context and reference files that are needed for that pipeline run, and dump them into the CRDS_PATH directory. Subsequent executions of the pipeline will first look to see if it has what it needs in CRDS_PATH and anything it doesn't have will be downloaded from the STScI cache. \n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-block alert-warning\">NOTE: Users at STScI should automatically have access to the Calibration Reference Data System (CRDS) cache for running the pipeline, and can skip setting these environment variables.</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Top of Notebook](#top)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='imports'></a>\n",
    "## Imports"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import packages necessary for this notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Packages that allow us to get information about objects:\n",
    "import asdf\n",
    "import copy\n",
    "import os\n",
    "import shutil\n",
    "\n",
    "# Numpy library:\n",
    "import numpy as np\n",
    "\n",
    "# For downloading data\n",
    "import requests\n",
    "\n",
    "# Astropy tools:\n",
    "from astropy.io import fits\n",
    "from astropy.utils.data import download_file\n",
    "from astropy.visualization import ImageNormalize, ManualInterval, LogStretch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set up matplotlib for plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "\n",
    "# Use this version for non-interactive plots (easier scrolling of the notebook)\n",
    "%matplotlib inline\n",
    "\n",
    "# Use this version (outside of Jupyter Lab) if you want interactive plots\n",
    "#%matplotlib notebook\n",
    "\n",
    "# These gymnastics are needed to make the sizes of the figures\n",
    "# be the same in both the inline and notebook versions\n",
    "%config InlineBackend.print_figure_kwargs = {'bbox_inches': None}\n",
    "\n",
    "mpl.rcParams['savefig.dpi'] = 80\n",
    "mpl.rcParams['figure.dpi'] = 80"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import JWST pipeline-related modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# List of possible data quality flags\n",
    "from jwst.datamodels import dqflags\n",
    "\n",
    "# The entire calwebb_detector1 pipeline\n",
    "from jwst.pipeline import calwebb_detector1\n",
    "\n",
    "# Individual steps that make up calwebb_detector1\n",
    "from jwst.dq_init import DQInitStep\n",
    "from jwst.saturation import SaturationStep\n",
    "from jwst.superbias import SuperBiasStep\n",
    "from jwst.ipc import IPCStep                                                                                    \n",
    "from jwst.refpix import RefPixStep                                                                \n",
    "from jwst.linearity import LinearityStep\n",
    "from jwst.persistence import PersistenceStep\n",
    "from jwst.dark_current import DarkCurrentStep\n",
    "from jwst.jump import JumpStep\n",
    "from jwst.ramp_fitting import RampFitStep\n",
    "from jwst import datamodels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check which version of the pipeline we are running:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import jwst\n",
    "print(jwst.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='convenience_functions'></a>\n",
    "## Define convenience functions and parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we define some functions and some parameters that we will use repeatedly throughout the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To make everything easier, all files saved by the pipeline\n",
    "# and pipeline steps will be saved to the working directory.\n",
    "# This will be more important for the level 2 and 3 pipelines\n",
    "# once we are working with association files.\n",
    "output_dir = './'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make sure the output directory exists before downloading any data\n",
    "if not os.path.exists(output_dir):\n",
    "    os.makedirs(output_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def download_files(files, output_directory):\n",
    "    \"\"\"Given a tuple or list of tuples containing (URL, filename),\n",
    "    download the given files into the current working directory.\n",
    "    Downloading is done via astropy's download_file. A symbolic link\n",
    "    is created in the specified output dirctory that points to the\n",
    "    downloaded file.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    files : tuple or list of tuples\n",
    "        Each 2-tuple should contain (URL, filename), where\n",
    "        URL is the URL from which to download the file, and\n",
    "        filename will be the name of the symlink pointing to\n",
    "        the downloaded file.\n",
    "        \n",
    "    output_directory : str\n",
    "        Name of the directory in which to create the symbolic\n",
    "        links to the downloaded files\n",
    "                \n",
    "    Returns\n",
    "    -------\n",
    "    filenames : list\n",
    "        List of filenames corresponding to the symbolic links\n",
    "        of the downloaded files\n",
    "    \"\"\"\n",
    "    # In the case of a single input tuple, make it a\n",
    "    # 1 element list, for consistency.\n",
    "    filenames = []\n",
    "    if isinstance(files, tuple):\n",
    "        files = [files]\n",
    "        \n",
    "    for file in files:\n",
    "        filenames.append(file[1])\n",
    "        if not os.path.isfile(os.path.join(output_directory, file[1])):\n",
    "            print('Downloading {}...'.format(file[1]))\n",
    "            demo_file = download_file(file[0], cache=True)\n",
    "            # Make a symbolic link using a local name for convenience\n",
    "            os.symlink(demo_file, os.path.join(output_directory, file[1]))\n",
    "        else:\n",
    "            print('{} already exists, skipping download...'.format(file[1]))\n",
    "            continue\n",
    "    return filenames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_jump(signal, jump_group, xpixel=None, ypixel=None, slope=None):\n",
    "    \"\"\"Function to plot the signal up the ramp for a\n",
    "    pixel and show the location of flagged jumps.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    signal : numpy.ndarray\n",
    "        1D array of signal values\n",
    "        \n",
    "    jump_group : list\n",
    "        List of boolean values whether a jump is present or\n",
    "        not in each group\n",
    "        \n",
    "    slope : numpy.ndarray\n",
    "        1D array of signal values constructed from the slope\n",
    "    \"\"\"\n",
    "    groups = np.arange(len(signal))\n",
    "    fig = plt.figure(figsize=(6, 6))\n",
    "    ax = plt.subplot()\n",
    "\n",
    "    plt.plot(groups, signal, marker='o', color='black')\n",
    "    plt.plot(groups[jump_group], signal[jump_group], marker='o', color='red',\n",
    "             label='Flagged Jump')\n",
    "    \n",
    "    if slope is not None:\n",
    "        plt.plot(groups, slope, marker='o', color='blue', label='Data from slope')\n",
    "        \n",
    "    plt.legend(loc=2)\n",
    "\n",
    "    plt.xlabel('Groups')\n",
    "    plt.ylabel('Signal (DN)')\n",
    "    fig.tight_layout()\n",
    "    plt.subplots_adjust(top=0.95)\n",
    "    \n",
    "    if xpixel and ypixel:\n",
    "        plt.title('Pixel ('+str(xpixel)+','+str(ypixel)+')')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_jumps(signals, jump_groups, pixel_loc, slopes=None):\n",
    "    \"\"\"Function to plot the ramp and show the jump location\n",
    "    for several pixels. For simplicity, let's force the input\n",
    "    number of pixels to be a square. \n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    signals : numpy.ndarray\n",
    "        2D array (groups x pix) of signal values\n",
    "        \n",
    "    jump_groups : numpy.ndarray\n",
    "        2D array containing boolean entries for each group of\n",
    "        each pixel, describing where the jumps were found\n",
    "        \n",
    "    pixel_loc : list\n",
    "        List of 2-tuples containing the (x, y)\n",
    "        location of the pixels with the jumps\n",
    "        \n",
    "    slopes : numpy.ndarray\n",
    "        2D array (groups x pix) of linear signal values\n",
    "        If not None, these will be overplotted onto the\n",
    "        plots of signals\n",
    "    \"\"\"\n",
    "    num_group, num_pix = signals.shape\n",
    "    root = np.sqrt(num_pix)\n",
    "    if int(root + 0.5) ** 2 != num_pix:\n",
    "        raise ValueError('Number of pixels input should be a square.')\n",
    "    \n",
    "    root = int(root)\n",
    "    groups = np.arange(num_group)\n",
    "    fig, axs = plt.subplots(root, root, figsize=(10, 10))\n",
    "\n",
    "    for index in range(len(pixel_loc)):\n",
    "        i = int(index % root)\n",
    "        j = int(index / root)\n",
    "        axs[i, j].plot(groups, signals[:, index], marker='o', color='black')\n",
    "        j_grp = jump_groups[:, index]\n",
    "        axs[i, j].plot(groups[j_grp], signals[j_grp, index],\n",
    "                       marker='o', color='red')\n",
    "        \n",
    "        if slopes is not None:\n",
    "            axs[i, j].plot(groups, slopes[:, index], marker='o', color='blue')\n",
    "        \n",
    "        axs[i, j].set_title('Pixel ({}, {})'.format(pixel_loc[index][1], pixel_loc[index][0]))\n",
    "        \n",
    "    plt.xlabel('Groups')\n",
    "    plt.ylabel('Signal (DN)')\n",
    "    fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_ramp(groups, signal, xpixel=None, ypixel=None, title=None):\n",
    "    \"\"\"Function to plot the up the ramp signal for a pixel.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    groups : numpy.ndarray\n",
    "        1D array of group numbers. X-axis values.\n",
    "        \n",
    "    signal : numpy.ndarray\n",
    "        1D array of pixel signal values.\n",
    "        \n",
    "    xpixel : int\n",
    "        X-coordinate of the pixel being plotted. Used for legend only.\n",
    "        \n",
    "    ypixel : int\n",
    "        Y-coordinate of the pixel being plotted. Used for legend only.\n",
    "        \n",
    "    title : str\n",
    "        String to use for the plot title\n",
    "    \"\"\"    \n",
    "    fig = plt.figure(figsize=(8, 8))\n",
    "    ax = plt.subplot()\n",
    "    if xpixel and ypixel:\n",
    "            plt.plot(groups, signal, marker='o',\n",
    "                     label='Pixel ('+str(xpixel)+','+str(ypixel)+')') \n",
    "            plt.legend(loc=2)\n",
    "\n",
    "    else:\n",
    "        plt.plot(groups, signal, marker='o')\n",
    "        \n",
    "    plt.xlabel('Groups')\n",
    "    plt.ylabel('Signal (DN)')\n",
    "    fig.tight_layout()\n",
    "    plt.subplots_adjust(left=0.15)\n",
    "    \n",
    "    if title:\n",
    "        plt.title(title)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_ramps(groups, signal1, signal2, label1=None, label2=None, title=None):\n",
    "    \"\"\"Function to plot the up the ramp signal for two pixels\n",
    "    on a single plot.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    groups : numpy.ndarray\n",
    "        1D array of group numbers. X-axis values.\n",
    "        \n",
    "    signal1 : numpy.ndarray\n",
    "        1D array of signal values for first pixel\n",
    "        \n",
    "    signal2 : numpy.ndarray\n",
    "        1D array of signal values for second pixel\n",
    "        \n",
    "    label1 : str\n",
    "        Label to place in the legend for pixel1\n",
    "        \n",
    "    label2 : str\n",
    "        Label to place in the legend for pixel2\n",
    "        \n",
    "    title : str\n",
    "        String to place in the title of the plot\n",
    "    \"\"\"    \n",
    "    fig = plt.figure(figsize=(6, 6))\n",
    "    ax = plt.subplot()\n",
    "    if label1:\n",
    "        plt.plot(groups, signal1, marker='o', color='black', label=label1)\n",
    "    else:\n",
    "        plt.plot(groups, signal1, marker='o', color='black')\n",
    "    if label2:\n",
    "        plt.plot(groups, signal2, marker='o', color='red', label=label2)\n",
    "    else:\n",
    "        plt.plot(groups, signal2, marker='o', color='red')\n",
    "    if label1 or label2:\n",
    "        plt.legend(loc=2)\n",
    "        \n",
    "    plt.xlabel('Groups')\n",
    "    plt.ylabel('Signal (DN)')\n",
    "    fig.tight_layout()\n",
    "    plt.subplots_adjust(left=0.15)\n",
    "    plt.subplots_adjust(top=0.95)\n",
    "    \n",
    "    if title:\n",
    "        plt.title(title)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def show_image(data_2d, vmin, vmax, xpixel=None, ypixel=None, title=None):\n",
    "    \"\"\"Function to generate a 2D, log-scaled image of the data, \n",
    "    with an option to highlight a specific pixel (with a red dot).\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    data_2d : numpy.ndarray\n",
    "        Image to be displayed\n",
    "        \n",
    "    vmin : float\n",
    "        Minimum signal value to use for scaling\n",
    "        \n",
    "    vmax : float\n",
    "        Maximum signal value to use for scaling\n",
    "        \n",
    "    xpixel : int\n",
    "        X-coordinate of pixel to highlight\n",
    "        \n",
    "    ypixel : int\n",
    "        Y-coordinate of pixel to highlight\n",
    "        \n",
    "    title : str\n",
    "        String to use for the plot title\n",
    "    \"\"\"\n",
    "    norm = ImageNormalize(data_2d, interval=ManualInterval(vmin=vmin, vmax=vmax),\n",
    "                          stretch=LogStretch())\n",
    "    fig = plt.figure(figsize=(8, 8))\n",
    "    ax = fig.add_subplot(1, 1, 1)\n",
    "    im = ax.imshow(data_2d, origin='lower', norm=norm)\n",
    "    \n",
    "    if xpixel and ypixel:\n",
    "        plt.plot(xpixel, ypixel, marker='o', color='red', label='Selected Pixel')\n",
    "\n",
    "    fig.colorbar(im, label='DN')\n",
    "    plt.xlabel('Pixel column')\n",
    "    plt.ylabel('Pixel row')\n",
    "    if title:\n",
    "        plt.title(title)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def side_by_side(data1, data2, vmin, vmax, title1=None, title2=None, title=None):\n",
    "    \"\"\"Show two images side by side for easy comparison. Optionally highlight\n",
    "    a given pixel with a red dot.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    data1 : numpy.ndarray\n",
    "        First image to be displayed\n",
    "        \n",
    "    data2 : numpy.ndarray\n",
    "        Second image to be displayed\n",
    "        \n",
    "    vmin : float\n",
    "        Minimum signal value to use for scaling\n",
    "        \n",
    "    vmax : float\n",
    "        Maximum signal value to use for scaling\n",
    "            \n",
    "    title1 : str\n",
    "        Title to use for first (left) plot\n",
    "        \n",
    "    title2 : str\n",
    "        Title to use for the second (right) plot\n",
    "\n",
    "    title : str\n",
    "        String to use for the plot title\n",
    "    \"\"\"\n",
    "    norm = ImageNormalize(data1, interval=ManualInterval(vmin=vmin, vmax=vmax),\n",
    "                          stretch=LogStretch())\n",
    "\n",
    "    fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(11, 8))\n",
    "    im = axes[0].imshow(data1, origin='lower', norm=norm)\n",
    "    im = axes[1].imshow(data2, origin='lower', norm=norm)\n",
    "    \n",
    "    axes[0].set_xlabel('Pixel column')\n",
    "    axes[0].set_ylabel('Pixel row')\n",
    "    axes[1].set_xlabel('Pixel column')\n",
    "    \n",
    "    if title1:\n",
    "        axes[0].set_title(title1)\n",
    "    if title2:\n",
    "        axes[1].set_title(title2)\n",
    "        \n",
    "    fig.subplots_adjust(right=0.8)\n",
    "    cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])\n",
    "    fig.colorbar(im, cax=cbar_ax, label='DN')\n",
    "    \n",
    "    if title:\n",
    "        fig.suptitle(title)    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Top of Notebook](#top)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='download_data'></a>\n",
    "## Download Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this module, we will use an uncalibrated NIRCam simulated imaging exposure that is stored in Box. Let's grab it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "uncal_info = ('https://stsci.box.com/shared/static/j46wpyirlbqo30e7c9719ycnuc1qk2lu.fits',\n",
    "              'jw98765001001_01101_00003_nrcb5_uncal.fits')\n",
    "uncal_file = download_files(uncal_info, output_dir)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Also download a \"trapsfilled\" file, which specifies the state of the charge traps at\n",
    "# the time of the preceding exposure. This will serve as input to the persistence\n",
    "# step.\n",
    "persist_info = ('https://stsci.box.com/shared/static/ehkof12d43h6nnpijs2r4tyuzde3nzc9.fits',\n",
    "                'jw98765001001_01101_00002_nrcb5_trapsfilled.fits')\n",
    "persist_file = download_files(persist_info, output_dir)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download example parameter reference files\n",
    "param_info = [('https://stsci.box.com/shared/static/nmaipgvlferrz95eyksra4zmlhx2m5bm.asdf',\n",
    "               'detector1pipeline_modified_paramfile.asdf'),\n",
    "              ('https://stsci.box.com/shared/static/kww0o8d77h2u03b38gtdr34g8vqygx8w.asdf',\n",
    "               'refpix_paramfile.asdf'),\n",
    "              ('https://stsci.box.com/shared/static/c95paaz9gsvrwwo3slx8k9qav242idi1.asdf',\n",
    "               'refpix_no_side_pix_paramfile.asdf')]\n",
    "detector1_param_reffile, refpix_param_reffile, refpix_param_reffile_no_side_pix = \\\n",
    "    download_files(param_info, output_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also download the example MIRI file if you wish to try the exercise of running it through the pipeline."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download raw MIRI file for notebook exercises\n",
    "miri_uncal_info = ('https://stsci.box.com/shared/static/rsav0kkg4dhb4y3oth6c0orcafc83qgm.fits',\n",
    "                   'miri_F770W_exp1_uncal.fits')\n",
    "miri_uncal_file = download_files(miri_uncal_info, output_dir)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download a MIRI parameter refrence file for the pipeline\n",
    "miri_param_info = ('https://stsci.box.com/shared/static/cbsj8v6ull992n4jdkzlg5485arsmars.asdf',\n",
    "                   'miri_detector1_param_file.asdf')\n",
    "miri_det1_param_reffile = download_files(miri_param_info, output_dir)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='calling_methods'></a>\n",
    "## Methods for calling steps/pipelines"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are three common methods by which the pipeline or pipeline steps can be called. From within python, the `run()` and `call()` methods of the pipeline or step classes can be used. Alternatively, the `strun` command can be used from the command line. Within this notebook, in the section where we [call the entire pipeline](#detector1_at_once), as well as the section where we [call the Reference Pixel Subtraction](#refpix) step, we show examples of all three methods. For the remainder of the pipeline steps, we will focus on using the `run()` method.\n",
    "\n",
    "When using the `call()` method or `strun`, optional input parameters can be specified via [parameter reference files](#parameter_reffiles). When using the `run()` method, these parameters are instead specified within python."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='parameter_reffiles'></a>\n",
    "## Parameter Reference Files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When calling a pipeline or pipeline step using the `call()` method or the command line, [parameter reference files](https://jwst-pipeline.readthedocs.io/en/stable/jwst/stpipe/config_asdf.html#config-asdf-files) can be used to specify values for input parameters. These reference files are in [asdf](https://asdf.readthedocs.io/en/stable/) format and appear somewhat similar to json files when examined in a text editor. \n",
    "\n",
    "Versions of parameter reference files containing default parameter values for each step and pipeline are available in CRDS. When using the `call()` method, if you do not specify a parameter reference file name in the call, the pipeline or step will retrieve and use the appropriate file from CRDS, which will then run the pipeline or step with the parameter values in that file. If you provide the name of a parameter reference file, then the parameter values in that file will take precedence. For any parameter not specified in your parameter reference file, the pipeline will use the default value.\n",
    "\n",
    "When using `strun`, the parameter reference file is a required input."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at the contents of a parameter reference file. We'll open it using the asdf package, and use the `tree` attribute to see what's inside:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "det1_reffile = asdf.open(detector1_param_reffile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "det1_reffile.tree"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The top part of the file contains various metadata entries about the file itself. Below that, you'll see a `'name'` entry, which lists `Detector1Pipeline` as the class to which these parameters apply. The next line contains the `parameters` entry, which lists parameters and values attached to the pipeline itself. Below this is the `steps` entry, which contains a list of dictionaries. Each dictionary refers to one step within the pipeline, and specifies parameters and values that apply to that step. If you look through these entries, you'll see the same parameters and values that we specified manually when using the `run()` method above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Don't forget to close the file\n",
    "det1_reffile.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Top of Notebook](#top)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "<a id='detector1'></a>\n",
    "## The calwebb_detector1 pipeline: Ramps to slopes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the sections below, we will run the Stage 1 pipeline on a single uncalibrated NIRCam file. We will first call the entire [*calwebb_detector1* pipeline](https://jwst-pipeline.readthedocs.io/en/latest/jwst/pipeline/calwebb_detector1.html#calwebb-detector1) on the file. The pipeline is a wrapper which will string together all of the appropriate steps in the proper order. \n",
    "\n",
    "The final output from this call is an uncalibrated slope image which is ready to go into the Stage 2 pipeline. \"Uncalibrated\" in this case means that the data are in units of DN/sec. In Stage 2 the flux calibration will be applied, at which point the data will be in physical units (e.g. MJy/str) and referred to as \"calibrated\".\n",
    "\n",
    "After that, we will go back to the original uncalibrated ramp and manually run it through each of the steps that comprise the Stage 1 pipeline. For each step we will describe in more detail what is going on, and for several of the more interesting steps, we will examine the output.\n",
    "\n",
    "See [Figure 1](https://jwst-docs.stsci.edu/jwst-data-reduction-pipeline/stages-of-processing/calwebb_detector1) on the calwebb_detector1 algorithm page for a map of which steps are performed on NIR data and which are used for MIRI data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_file_base = os.path.basename(uncal_file).replace('uncal.fits', '')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='detector1_at_once'></a>\n",
    "### Run the entire `calwebb_detecor1` pipeline\n",
    "\n",
    "In this section we show how to run the entire calwebb_detector1 pipeline with a single call. In this case the pipeline code can determine which instrument was used to collect the data and runs the appropriate steps in the proper order.\n",
    "\n",
    "We set parameter values for some of the individual steps, save some outputs, etc, and then call the pipeline.\n",
    "\n",
    "We will call the pipeline using the `run()` method and while that is running, we will go over the equivalent `call()` and `strun` commands, examine some of the pipeline log entries that are printed to the screen, and then look at the pipeline output."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='run_method'></a>\n",
    "#### Call the pipeline using the run() method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When using the `run()` method to execute a pipeline (or step), the pipeline class is first instantiated without the data to be processed. Optional input parameters are specified using attributes of the class instance. Finally, the call to the `run()` method is made and the data are supplied.  See here for [more examples of the run() method](https://jwst-pipeline.readthedocs.io/en/stable/jwst/stpipe/call_via_run.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `run()` method does not take any kind of parameter reference file as input. If you wish to set values for various parameters, you must do that manually. Below, we set several parameters in order to show how it's done. \n",
    "\n",
    "How do you know what parameters are available to be set and what their default values are? The `spec` property for individual steps will list them. The property is less useful for the pipelines themselves, as it does not show the parameters for the steps comprising the pipeline.\n",
    "\n",
    "All steps and pipelines have several common parameters that can be set. \n",
    "\n",
    "* `save_results` specifies whether or not to save the output of that step/pipeline to a file. The default is False.\n",
    "* `output_dir` is the directory into which the output files will be saved.\n",
    "* `output_file` is the base filename to use for the saved result. Note that each step/pipeline will add a custom suffix onto output_file. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at the available parameters for the reference pixel subtraction step, and the cosmic ray flagging step, and manually set some of these in our call to `run()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(RefPixStep.spec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(JumpStep.spec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='detector1_using_run'></a>\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "Finally, let's run the pipeline. The output can be a little overwhelming. There will be multiple log entries printed to the screen for each step.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Create an instance of the pipeline class\n",
    "detector1 = calwebb_detector1.Detector1Pipeline()\n",
    "\n",
    "# Set some parameters that pertain to the\n",
    "# entire pipeline\n",
    "detector1.output_dir = output_dir\n",
    "detector1.save_results = True\n",
    "\n",
    "# Set some parameters that pertain to some of\n",
    "# the individual steps\n",
    "detector1.refpix.use_side_ref_pix = True\n",
    "detector1.linearity.save_results = True\n",
    "detector1.jump.rejection_threshold = 6\n",
    "\n",
    "# Specify the name of the trapsfilled file, which\n",
    "# contains the state of the charge traps at the end\n",
    "# of the preceding exposure\n",
    "detector1.persistence.input_trapsfilled = persist_file\n",
    "\n",
    "# Call the run() method\n",
    "run_output = detector1.run(uncal_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This cell will take a few minutes to run. While we're waiting, let's look at the other two methods that can be used to call the pipeline. Then we'll come back here and look at the log meessages output by this cell so we can see what happened."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='call_method'></a>\n",
    "#### Call the pipeline using the call() method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When using the `call()` method, a single command will instantiate and run the pipeline (or step). The input data and optional parameter reference files are supplied in this single command. In this case, any desired input parameters cannot be set after instantiation, as with the `run()` method. See here for [example usage of call() method](https://jwst-pipeline.readthedocs.io/en/stable/jwst/stpipe/call_via_call.html).\n",
    "\n",
    "The commands below will call the pipeline using the `call()` method and will supply the parameter reference file. Since we just ran the pipeline with the `run()` method above, we won't actually execute the call to `call()`. But if you wish to try it out, use the pull-down menu above to change the cell to be 'Code', and then execute it. (Or, Click 'Cell' > 'Cell Type' > 'Code')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "\n",
    "<b>Method #1:</b>\n",
    "Provide the name of the observation file, the pipeline-specific input paramters, and the name of the parameter reference file that specifies step-specific parameters\n",
    "</div>"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# NOTE: This is a \"raw\" cell and will not run unless you change to \"code\"\n",
    "call_output = calwebb_detector1.Detector1Pipeline.call(uncal_file,\n",
    "                                                      output_dir=output_dir,\n",
    "                                                      save_results=True,\n",
    "                                                      config_file=detector1_param_reffile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "\n",
    "<b>Method #2:</b>\n",
    "In this case, build a nested dictionary that specifies parameter values for various steps, and provide it in the call to call().\n",
    "</div>"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# NOTE: This is a \"raw\" cell and will not run unless you change to \"code\"\n",
    "parameter_dict = {\"refpix\": {\"use_side_ref_pixels\": True},\n",
    "                  \"linearity\": {\"save_results\": True},\n",
    "                  \"jump\": {\"rejection_threshold\": 6},\n",
    "                 }\n",
    "call_output = calwebb_detector1.Detector1Pipeline.call(uncal_file, output_dir=output_dir, save_results=True,\n",
    "                                                       steps=parameter_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='command_line'></a>\n",
    "#### Call the pipeline using the command line"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calling a pipeline or step from the command line is similar to using the `call()` method. The data file to be processed, along with an optional parameter reference file and optional parameter/value pairs can be provided to the `strun` command. See here for [additional examples of command line calls](https://jwst-pipeline.readthedocs.io/en/stable/jwst/introduction.html?highlight=%22command%20line%22#running-from-the-command-line).\n",
    "\n",
    "The cells below contains two different command line commands that use `strun` to call the calwebb_detector1 pipeline. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "\n",
    "<b>Method #1:</b>\n",
    "We provide the name of the pipeline class, the observation file, and explicitly set pipeline- and step-specific parameters. You can see that the command quickly becomes quite large with the added parameter settings. \n",
    "    \n",
    "```\n",
    "    strun jwst.pipeline.Detector1Pipeline jw98765001001_01101_00003_nrcb5_uncal.fits --save_results=True --output_dir='./' --steps.refpix.use_side_ref_pixels=True --steps.linearity.save_results=True --steps.jump.rejection_threshold=6 \n",
    "```\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "\n",
    "<b>Method #2:</b>\n",
    "This version of the command is much more succinct, as the parameter values to be set are all contained within the parameter reference file. The pipeline class is also contained in the parameter reference file, so there is no need to specify it in the command itself.\n",
    "    \n",
    "```\n",
    "    strun detector1pipeline_modified_paramfile.asdf jw98765001001_01101_00003_nrcb5_uncal.fits\n",
    "```\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Examine the outputs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The primary output of the calwebb_detector1 pipeline is a file containing a rate image for the exposure. The units of the data are ADU/sec.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rate_file = uncal_file.replace('uncal.fits', 'rate.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rate_file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rate_data = fits.getdata(rate_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_image(rate_data, 0.5, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also, since we set the linearity step's `save_results` parameter to True in our calls above, the pipeline saved the output of the linearity step. In this case, the output file will have the same name as the input uncal file, but with the suffix 'linearity' rather than 'uncal'. \n",
    "\n",
    "**NOTE:** This differs slightly from the case where we call the linearity step itself and save the results. In that case, as we will see in the [linearity](#linearity) section, the output file will have the suffix 'linearitystep' rather than 'linearity'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "linear_file = rate_file.replace('rate.fits', 'linearity.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lin_data = fits.getdata(linear_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will look in more detail at the effects of the linearity correction step in the [linearity](#linearity) section below. For now, let's just look at the final group of the integration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's look at the data in the final group of the linearized data:\n",
    "show_image(lin_data[0, -1, :, :], 100, 10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Top of Notebook](#top)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='exercise1'></a>\n",
    "### Exercise 1: Run calwebb_detector1 on MIRI data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try running the pipeline on the raw MIRI file downloaded earlier. We downloaded a parameter reference file along with it (**miri_detector1_param_file.asdf**), so you can try either the `run()` or `call()` methods. Note that MIRI data go through some additional steps compared to those we saw for the NIRCam data above. There aren't many parameters beyond those from the NIRCam example to adjust, but feel free to try tweaking some values. Also, try saving the output from the dark current subtraction step."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See [Exercise 1 solutions](#exercise1_solution)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using the run() method:\n",
    "# Instantiate the pipeline\n",
    "\n",
    "# Save the final output of the pipeline\n",
    "\n",
    "# Save the output from the dark current subtraction step,\n",
    "# and let's use a more stringent limit for cosmic ray\n",
    "# detection\n",
    "\n",
    "# Call the run() method\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using the call() method:\n",
    "# First edit the parameter reference file:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the call() method and supply the parameter reference file\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Look at one of the output rate files\n",
    "# HINT: Open one of the rate files and use show_image with\n",
    "# min and max signal rate limits of 0 and 100.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='detector1_step_by_step'></a>\n",
    "## Run the individual pipeline steps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the sections below we run the steps contained within calwebb_detector1 one at a time, in order to more clearly see what each step is doing. Since our example data are from NIRCam, we will skip the MIRI-specific steps."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='dq_init'></a>\n",
    "### The `Data Quality Initialization` step\n",
    "\n",
    "#### Summary\n",
    "\n",
    "This step populates the Data Quality (DQ) mask that is associated with the data file. The DQ flags from the `MASK` reference file are copied into the `PIXELDQ` extension of the input file. A table showing the [mapping of bit values](https://jwst-pipeline.readthedocs.io/en/stable/jwst/references_general/references_general.html#data-quality-flags) in the `MASK` file decribes what types of bad pixels can be flagged. Any other bad pixel types will be ignored.\n",
    "\n",
    "#### Documentation\n",
    "\n",
    "[Full description](https://jwst-pipeline.readthedocs.io/en/stable/jwst/dq_init/description.html) of the step.\n",
    "\n",
    "#### Arguments\n",
    "\n",
    "There are no optional arguments for this step\n",
    "\n",
    "#### Reference files used\n",
    "\n",
    "This step uses the `MASK` reference file. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Run the step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Using the run() method\n",
    "dq_init_step = DQInitStep()\n",
    "dq_init_step.output_dir = output_dir\n",
    "dq_init_step.save_results = True\n",
    "\n",
    "# Note that the run() method can be called on EITHER:\n",
    "# the datamodel instance output from the group_scale\n",
    "# step.\n",
    "dq_init = dq_init_step.run(uncal_file)\n",
    "\n",
    "# OR:\n",
    "# the fits file containing the group_scale output\n",
    "#dq_init = dq_init_step.run(group_scale_output_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The step finished without crashing, but as it is said above, there are some warnings worth noting.\n",
    "Notably, the `NOISY` and `WEIRD` do not correspond to existing `DQ` mnemonics, so they are ignored. This is expected, and means that the `MASK` reference file contains some pixels flagged as `NOISY` and `WEIRD`. Since these bad pixel types are not present in the list of known types of bad pixels, as shown below, these flags are ignored."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dqflags.pixel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The pixel values in the `SCI` extension are not changed in this step. Instead, the DQ flags are copied into the `PIXELDQ` extension. The `GROUPDQ` values are not changed in this step. Let's check the `PIXELDQ` values and see what has changed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have the datamodel instance of the output available already, but if you wanted to open the output file from the data quality initiailzation step, the output file has the same name as the input file, with the original suffix replaced by \"dqinitstep\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dq_init_output_file = os.path.join(output_dir, '{}dqinitstep.fits'.format(input_file_base))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx_pixelDQ = np.where(dq_init.pixeldq.flatten() == 0.)[0]\n",
    "num_flagged = dq_init.pixeldq.size - len(idx_pixelDQ)\n",
    "print('Total pixels in PIXELDQ: {}'.format(dq_init.pixeldq.size))\n",
    "print('{} pixels have no flags.'.format(len(idx_pixelDQ)))\n",
    "print('{} pixels ({:.2f}% of the detector) have flags.'.format(num_flagged, num_flagged / dq_init.pixeldq.size))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Top of Notebook](#top)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='saturation'></a>\n",
    "## The `Saturation Flagging` step\n",
    "\n",
    "#### Summary\n",
    "\n",
    "This step checks the signal values in all pixels across all groups, and adds a [`saturated` flag](https://jwst-pipeline.readthedocs.io/en/stable/jwst/references_general/references_general.html#data-quality-flags) to the `GROUPDQ` extension for pixels and groups where the signal is above the saturation limit.\n",
    "\n",
    "#### Documentation\n",
    "\n",
    "[Full description](https://jwst-pipeline.readthedocs.io/en/stable/jwst/saturation/description.html) of the step.\n",
    "\n",
    "#### Arguments\n",
    "\n",
    "There are no optional arguments for this step\n",
    "\n",
    "#### Reference files used\n",
    "\n",
    "This step uses the [`SATURATION`](https://jwst-pipeline.readthedocs.io/en/stable/jwst/saturation/reference_files.html) reference file. This file contains a map of the saturation threshold in ADU for each pixel on the detector."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Run the step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Using the run() method\n",
    "saturation_step = SaturationStep()\n",
    "saturation_step.output_dir = output_dir\n",
    "saturation_step.save_results = True\n",
    "\n",
    "# Call using the the output from the previously-run dq_init step\n",
    "saturation = saturation_step.run(dq_init)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If there are any saturated values, they should appear in the `GROUPDQ` arrays. Let's examine the `GROUPDQ` data and see if there are any detected:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "saturated = np.where(saturation.groupdq & dqflags.pixel['SATURATED'] > 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_sat_flags = len(saturated[0])\n",
    "print(('Found {} saturated flags. This may include multiple saturated '\n",
    "       'groups within a given pixel'.format(num_sat_flags)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's find a pixel that saturated part of the way up the ramp."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a 4D boolean map of whether the saturation flag is present or not.\n",
    "saturated = (saturation.groupdq & dqflags.pixel['SATURATED'] > 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Collapse that down to a 2D map that lists the number of saturated groups \n",
    "# for each pixel.\n",
    "saturated_2d = np.sum(saturated[0, :, :, :], axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get coordinates of pixels that are saturated in some, but not all, groups.\n",
    "partial_sat = np.where((saturated_2d > 0) & (saturated_2d < saturated.shape[1]))\n",
    "print(\"{} pixels are partially saturated.\".format(len(partial_sat[0])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's choose one of these partially saturated pixels and look at the signal values up the ramp, along with which groups are flagged as saturated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sat_y, sat_x = partial_sat\n",
    "sat_index = 123\n",
    "y = sat_y[sat_index]\n",
    "x = sat_x[sat_index]\n",
    "grps = saturated[0, :, y, x]\n",
    "print('Saturation flags up the ramp (0 is not saturated, 2 is saturated): {}'\n",
    "      .format(saturation.groupdq[0, :, y, x]))\n",
    "print('Pixel signal values up the ramp: {}'.format(saturation.data[0, :, y, x]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot these in order to get a clearer look at the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "groups = np.arange(saturation.data.shape[1])\n",
    "full_ramp = saturation.data[0, :, y, x]\n",
    "sat_dq = saturation.groupdq[0, :, y, x].astype(bool)\n",
    "saturated_points = copy.deepcopy(saturation.data[0, :, y, x])\n",
    "saturated_points[~sat_dq] = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_ramps(groups, full_ramp, saturated_points, label1='Not Saturated',\n",
    "           label2='Saturated', title='Pixel ({}, {})'.format(x, y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " <a id='superbias'> </a>\n",
    "## The `Superbias Subtraction` step\n",
    "\n",
    "#### Summary\n",
    "\n",
    "This step subtracts the superbias reference frame from each group of the science exposure.\n",
    "\n",
    "#### Documentation\n",
    "\n",
    "[Full description](https://jwst-pipeline.readthedocs.io/en/stable/jwst/superbias/description.html) of the step.\n",
    "\n",
    "#### Arguments\n",
    "\n",
    "There are no optional arguments for this step\n",
    "\n",
    "#### Reference files used\n",
    "\n",
    "This step uses the [`SUPERBIAS`](https://jwst-pipeline.readthedocs.io/en/stable/jwst/superbias/reference_files.html) reference file. This file contains a map of the superbias signal in ADU for each pixel on the detector.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Run the step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Using the run() method\n",
    "superbias_step = SuperBiasStep()\n",
    "superbias_step.output_dir = output_dir\n",
    "superbias_step.save_results = True\n",
    "\n",
    "# Call using the the output from the previously-run saturation step\n",
    "superbias = superbias_step.run(saturation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visually compare the science products to the raw `uncal` data, looking only at the last group of the first integration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "superbias_output_file = os.path.join(output_dir, '{}superbiasstep.fits'.format(input_file_base))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "superbias_output_file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "superbias.data[0, 0, :, :].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "side_by_side(saturation.data[0, 0, :, :], superbias.data[0, 0, :, :], vmin=10000, vmax=18000,\n",
    "            title1='Before superbias subtraction', title2='After superbias subtraction')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Top of Notebook](#top)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='refpix'></a>\n",
    "## The `Reference Pixel Subtraction` step\n",
    "\n",
    "#### Summary\n",
    "\n",
    "This step uses the reference pixels, which are not sensitive to illumination, to subtract group- and amplifier-dependent signal originating in the readout electronics from the data. There are two distinct corrections that are applied here.\n",
    "\n",
    "First, the rows of reference pixels on the top and bottom of the detector are used to subtract amplifier-dependent offsets from each group. Within a given amplifier, the mean value of all reference pixels in even numbered columns is subtracted from the science pixels in the even numbered columns. The same strategy is used for the odd numbered columns. \n",
    "\n",
    "The second part of the reference pixel subtraction step uses the reference pixels along the left and right sides of the detector to mitigate 1/f noise. This noise is visible in the data as horizontal banding that stretches across the entire detector. \n",
    "\n",
    "#### Documentation\n",
    "\n",
    "[Full description](https://jwst-pipeline.readthedocs.io/en/stable/jwst/refpix/description.html) of the step.\n",
    "\n",
    "#### Arguments\n",
    "\n",
    "[Full details on the optional arguments](https://jwst-pipeline.readthedocs.io/en/stable/jwst/refpix/arguments.html).\n",
    "\n",
    "#### Reference files used\n",
    "\n",
    "This step does not use any reference files.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Run the step"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to better show the effects from the 2 parts of the reference file subtraction, we're going to run this step twice. First we'll perform only the mean value subtraction using the top and bottom reference pixels. Then on the second run, we'll perform both the mean value subtraction and the 1/f subtraction."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='refpix_run'></a>\n",
    "##### Using the run() method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A reminder of the available parameters that can be set:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(RefPixStep.spec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For this 'partial' run, we need to turn off the 1/f correction that\n",
    "# uses the reference pixels on the sides of the detector. Also, we'll\n",
    "# save the output using a unique name so as not to confuse the file\n",
    "# with the output where we run the entire repix subtraction step.\n",
    "\n",
    "refpix_step_no_sidepix = RefPixStep()\n",
    "refpix_step_no_sidepix.output_dir = output_dir\n",
    "refpix_step_no_sidepix.save_results = True\n",
    "refpix_step_no_sidepix.output_file = 'refpix_test_no_side_pixels'\n",
    "\n",
    "# Turn off the 1/f correction\n",
    "refpix_step_no_sidepix.use_side_ref_pixels = False\n",
    "\n",
    "# Call using the the output from the previously-run superbias step\n",
    "refpix_no_sidepix = refpix_step_no_sidepix.run(superbias)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next run the full correction. This will produce the output that we will feed into subsequent steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Instantiate and set parameters\n",
    "refpix_step = RefPixStep()\n",
    "refpix_step.output_dir = output_dir\n",
    "refpix_step.save_results = True\n",
    "\n",
    "# Call using the saturation instance from the previously-run\n",
    "# saturation step\n",
    "refpix = refpix_step.run(superbias)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='refpix_call'></a>\n",
    "##### Using the call() method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cells show how to run the reference pixel subtraction step with the same parameters as the preceeding two cells, but using the `call()` method rather than the `run()` method. In this case we have a separate parameter reference file for each call. \n",
    "\n",
    "Let's look at the parameter referece files. The difference between the two is in the `use_side_ref_pixels` entry, on the bottom line of each. These files look fairly similar to that for the Detector1Pipeline. In this case, the step's parameters and values are all listed in the `parameters` entry of the file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Open the parameter reference file and look at the tree\n",
    "default_refpix_params = asdf.open(refpix_param_reffile)\n",
    "default_refpix_params.tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Close the file\n",
    "default_refpix_params.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Open the parameter reference file where side reference pixels were not used\n",
    "# and look at the tree\n",
    "default_refpix_params = asdf.open(refpix_param_reffile_no_side_pix)\n",
    "default_refpix_params.tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Close the file\n",
    "default_refpix_params.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now use the `call()` method and the parameter reference file with all of the default settings to run the refernce pixel subtraction step. If you wish to run these cells, change the cell type to 'Code' in the pull-down menu above, and then execute."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# NOTE: This is a \"raw\" cell and will not run unless you change to \"code\"\n",
    "refpix_output = RefPixStep.call(superbias, output_dir=output_dir,\n",
    "                                save_results=True, \n",
    "                                config_file=refpix_param_reffile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Call the reference pixel subtraction step again, but use the parameter reference file that has the side reference pixel step turned off."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# NOTE: This is a \"raw\" cell and will not run unless you change to \"code\"\n",
    "refpix_output_no_side_pix = RefPixStep.call(superbias, output_dir=output_dir,\n",
    "                                            output_file='refpix_test_no_side_pixels_via_call',\n",
    "                                            save_results=True, \n",
    "                                            config_file=refpix_param_reffile_no_side_pix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='refpix_strun'></a>\n",
    "##### From the command line"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we show the commands that can be used from the command line to run the reference pixel correction step in the same ways as above. In this case, we would have had to save the output from the superbias subtraction step into a fits file, and then provide that file name in the command. For the purposes of this example, assume that the superbias subtraction output file is called superbias_sub.fits."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "    \n",
    "```\n",
    "strun refpix_paramfile.asdf superbias_sub.fits\n",
    "````\n",
    "    \n",
    "```\n",
    "strun refpix_no_side_pix_paramfile.asdf superbias_sub.fits\n",
    "```\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Examine the output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you wish to examine the output file from this step, the file is saved with the \"refpixstep\" suffix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the reference pixel step output filename\n",
    "refpix_output_file = os.path.join(output_dir, '{}refpixstep.fits'.format(input_file_base))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's explore the output from this step. As with all NIR detectors, the outermost 4 rows and columns comprise the reference pixels.\n",
    "\n",
    "Let's use the datamodels from before and after the reference pixel subtraction to have a look at the data. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll zoom in on the top few rows in order to see the changes. Note how the \"odd/even\" effect dominates the signal prior to reference pixel subtraction. Even numbered columns and odd numbered columns have significantly different signal levels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Side by side, before/after reference pixel subtraction\n",
    "side_by_side(superbias.data[0, 5, 2030:, 1030:1050], refpix.data[0, 5, 2030:, 1030:1050],\n",
    "             vmin=0, vmax=30000, title1='Before Refpix Subtraction',\n",
    "             title2='After Refpix Subtraction')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Image of the difference before/after reference pixel subtraction\n",
    "show_image(superbias.data[0, 5, 2030:, 1030:1050] - refpix.data[0, 5, 2030:, 1030:1050],\n",
    "           vmin=7000, vmax=16000, title='Difference: Before - After Refpix Subtraction')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's look at the difference in the data after the mean value subtraction compared to the case where both the mean value subtraction and the 1/f correction are done."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Side by side view of the reference pixel subtraction with and without\n",
    "# using the side reference pixels to subtract 1/f noise\n",
    "side_by_side(refpix_no_sidepix.data[0, 5, :, :], refpix.data[0, 5, :, :],\n",
    "             vmin=20, vmax=150, title1='No Side Repix Correction',\n",
    "             title2='With Side Refpix Correction')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Look at the difference image before and after the side reference\n",
    "# pixels are used to subtract the 1/f noise\n",
    "show_image(refpix.data[0, 5, :, :] - refpix_no_sidepix.data[0, 5, :, :],\n",
    "           vmin=-10, vmax=10, title=\"Difference with/without using side refpix\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Top of Notebook](#top)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='exercise2'></a>\n",
    "### Exercise 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try running the reference pixel subtraction step yourself! Try adjusting the number of rows of reference pixels that are smoothed as the step attempts to subtract 1/f noise. \n",
    "\n",
    "**HINT:** remember to use `RefPixStep.spec` to see the names of available parameters as well as their default values. Since we didn't explicitly set the smoothing length in our calls above, the default value was used. You can run the step using the `superbias` data model object, as we did above, or you can run on the `superbias_output_file` that we saved when we ran the SuperBiasStep."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See [Exercise 2 solution](#exercise2_solution)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check the names of the available parameters\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using the run() method.\n",
    "\n",
    "# Instantiate the step\n",
    "\n",
    "\n",
    "# Set the side_smoothing_length to a non-default value\n",
    "# (default is 11)\n",
    "\n",
    "\n",
    "# Call using the superbias instance from the previously-run\n",
    "# saturation step\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using the call() method\n",
    "\n",
    "# First open 'refpix_paramfile.asdf' and set the \n",
    "# side_smoothing_length entry to a new value.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='linearity'></a>\n",
    "## The `Linearity Correction` step \n",
    "\n",
    "#### Summary\n",
    "\n",
    "This step applies the classical linearity correction to the data on a pixel-by-pixel, integration-by-integration, group-by-group manner.\n",
    "\n",
    "#### Documentation\n",
    "\n",
    "[Full description](https://jwst-pipeline.readthedocs.io/en/stable/jwst/linearity/description.html) of the step.\n",
    "\n",
    "#### Arguments\n",
    "\n",
    "There are no optional arguments for this step\n",
    "\n",
    "#### Reference files used\n",
    "\n",
    "This step uses the [`LINEARITY`](https://jwst-pipeline.readthedocs.io/en/stable/jwst/linearity/reference_files.html) reference file. This file contains the polynomial coefficients used to apply the linearity correction to non-linear data.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Run the step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using the run() method\n",
    "linearity_step = LinearityStep()\n",
    "linearity_step.output_dir = output_dir\n",
    "linearity_step.save_results = True\n",
    "\n",
    "# Call using the refpix instance from the previously-run\n",
    "# refpix step\n",
    "linearity = linearity_step.run(refpix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Examine the output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output file from this step has the \"linearitystep\" suffix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "linearity_output_file = os.path.join(output_dir, '{}linearitystep.fits'.format(input_file_base))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at the signal up the ramp for a high signal pixel, in order to more easily see how the linearity correction changed the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using the 3rd group, find the difference between the data before and after\n",
    "# linearity correction. Find the pixels where this difference is greater than\n",
    "# 20 DN, and also where the signal in the final group is over 40,000 DN.\n",
    "lin_fix = linearity.data[0, 3, :, :] - refpix.data[0, 3, :, :]\n",
    "well_exposed = np.where((linearity.data[0, -1, :, :] > 40000.) & (lin_fix > 20))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('{} pixels meet the criteria above.'.format(len(well_exposed[0])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pick one of these pixels and plot the signal before and after the linearity correction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "index = 3\n",
    "lin_pix_x, lin_pix_y = (well_exposed[1][index], well_exposed[0][index])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "group_nums = np.arange(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_ramps(group_nums, refpix.data[0, :, lin_pix_y, lin_pix_x],\n",
    "           linearity.data[0, :, lin_pix_y, lin_pix_x], label1='Uncorrected', label2='Corrected',\n",
    "           title='Pixel ({}, {})'.format(lin_pix_x, lin_pix_y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the pixel reached saturation in group 9. So in group 9, the linearity correction made no changes. Between groups 0 to 8, you can see the original signal (in black) becoming more and more non-linear as signal increases, along with how the linearity correction modified the signal (red)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Top of Notebook](#top)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='persistence'></a>\n",
    "## The `Persistence Correction` step \n",
    "\n",
    "#### Summary\n",
    "\n",
    "This step uses a model to calculate the amount of signal in each group of each pixel that comes from persistence. This persistence signal is then subtracted, pixel-by-pixel and group-by-group, from the data.\n",
    "\n",
    "#### Documentation\n",
    "\n",
    "[Full description](https://jwst-pipeline.readthedocs.io/en/stable/jwst/persistence/description.html) of the step.\n",
    "\n",
    "#### Arguments\n",
    "\n",
    "[Optional arguments](https://jwst-pipeline.readthedocs.io/en/stable/jwst/persistence/arguments.html) for this step include setting the threshold signal value above which pixels are flagged in the DQ extension, as well as saving the subtracted persistence signal in a separate file.\n",
    "\n",
    "#### Reference files used\n",
    "\n",
    "This step uses the [`TRAPDENSITY`, `PERSAT`, and `TRAPPARS`](https://jwst-pipeline.readthedocs.io/en/stable/jwst/persistence/reference_files.html) reference files. The TRAPDENSITY file contains a map of the relative number of traps per pixel. The PERSAT reference file contains a map of the persistence saturation level, and the TRAPPARS reference file contains parameters related to the persistence calculation model. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Run the step"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "List the available parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(PersistenceStep.spec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using the run() method\n",
    "persist_step = PersistenceStep()\n",
    "persist_step.output_dir = output_dir\n",
    "persist_step.save_results = True\n",
    "\n",
    "# Specify the trapsfilled file, which contains the\n",
    "# state of the charge traps in the preceding exposure\n",
    "persist_step.input_trapsfilled = persist_file\n",
    "\n",
    "# Let's also save a separate file that contains the\n",
    "# subtracted persistence signal \n",
    "persist_step.save_persistence = True\n",
    "\n",
    "# Call using the refpix instance from the previously-run\n",
    "# linearity step\n",
    "persist = persist_step.run(linearity)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Examine the output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The primary output file from this step has the \"persistencestep.fits\" suffix added."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "persist_output_file = os.path.join(output_dir, '{}persistencestep.fits'.format(input_file_base))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at the optional output file that contains the subtracted persistence signal. \n",
    "\n",
    "<div class=\"alert alert-block alert-warning\">NOTE: In this case the output is pretty boring. It turns out that the trap density map reference file for NIRCam is empty because it is not yet well characterized. This means that the persistence signal is zero in all pixels for all exposures. Once the map is updated in commissioning, this step will start calculating persistence values.</dev>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "persist_signal_file = os.path.join(output_dir, '{}output_pers.fits'.format(input_file_base))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "persist_signal = fits.getdata(persist_signal_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The file contains the persistence signal for each group of each integration, as we can see by the shape."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "persist_signal.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.min(persist_signal), np.max(persist_signal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at the persistence signal in the final group."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_image(persist_signal[0, -1, :, :], 0, .01)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='dc'></a>\n",
    "##  The `Dark Current Subtraction` step \n",
    "\n",
    "#### Summary\n",
    "\n",
    "This step subtracts the dark current, group by group, from the input integrations.\n",
    "\n",
    "#### Documentation\n",
    "\n",
    "[Full description](https://jwst-pipeline.readthedocs.io/en/stable/jwst/dark_current/description.html) of the step.\n",
    "\n",
    "#### Arguments\n",
    "\n",
    "There are no optional arguments for this step\n",
    "\n",
    "#### Reference files used\n",
    "\n",
    "This step uses the [`DARK`](https://jwst-pipeline.readthedocs.io/en/stable/jwst/dark_current/reference_files.html) reference file. This file contains the measured mean dark current associated with the detector and subarray.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Run the step"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NOTE: This is one of the longer-running steps of calwebb_detector1. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using the run() method\n",
    "dark_step = DarkCurrentStep()\n",
    "dark_step.output_dir = output_dir\n",
    "dark_step.save_results = True\n",
    "\n",
    "# Call using the persistence instance from the previously-run\n",
    "# persistence step\n",
    "dark = dark_step.run(persist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dark_output_file = os.path.join(output_dir, '{}darkcurrentstep.fits'.format(input_file_base))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Top of Notebook](#top)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='jump'></a>\n",
    "## The `Cosmic Ray Flagging` step\n",
    "\n",
    "#### Summary\n",
    "\n",
    "This step searches for \"jumps\" in the ramp data. In this case, a jump in a pixel's ramp is defined as a large deviation in the count rate relative to that in the other groups. When a jump is found, the associated flag is added to the `GROUPDQ` extension for the group and pixel where the jump was detected. The science data are not modified at all. In the subsequent ramp-fitting step, the algorithm will look into the `GROUPDQ` array and ignore any groups where the jump flag has been set.\n",
    "\n",
    "#### Documentation\n",
    "\n",
    "[Full description](https://jwst-pipeline.readthedocs.io/en/stable/jwst/jump/description.html) of the step.\n",
    "\n",
    "#### Arguments\n",
    "\n",
    "The jump step has [several optional arguments](https://jwst-pipeline.readthedocs.io/en/stable/jwst/jump/arguments.html)\n",
    "\n",
    "#### Reference files used\n",
    "\n",
    "This step uses the [`READNOISE`](https://jwst-pipeline.readthedocs.io/en/stable/jwst/references_general/readnoise_reffile.html) and [`GAIN`](https://jwst-pipeline.readthedocs.io/en/stable/jwst/references_general/gain_reffile.html) reference files. These files contain maps of the readnoise and gain values across the detector.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Run the step"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "List the available parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(JumpStep.spec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using the run() method\n",
    "jump_step = JumpStep()\n",
    "jump_step.output_dir = output_dir\n",
    "jump_step.save_results = True\n",
    "jump_step.rejection_threshold = 9\n",
    "\n",
    "# Call using the dark instance from the previously-run\n",
    "# dark current subtraction step\n",
    "jump = jump_step.run(dark)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Examine the output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "jump_output_file = os.path.join(output_dir, '{}jumpstep.fits'.format(input_file_base))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see what some of these jumps look like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "jump_flags = np.where(jump.groupdq & dqflags.pixel['JUMP_DET'] > 0)\n",
    "print('{} jump flags detected.'.format(len(jump_flags[0])))\n",
    "\n",
    "jump_map = (jump.groupdq & dqflags.pixel['JUMP_DET'] > 0)\n",
    "jump_map_2d = np.sum(jump_map[0, :, :, :], axis=0)\n",
    "jump_map_indexes = np.where(jump_map_2d > 0)\n",
    "impacted_pix = np.sum(jump_map_2d > 0)\n",
    "total_pix = 2048 * 2048\n",
    "print(('{} pixels ({:.2f}% of the detector) have been flagged with '\n",
    "      'at least one jump.'.format(impacted_pix, 100. * impacted_pix / total_pix)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot some of the pixels impacted by jumps. The red marks signify flagged jumps. These signal values will be ignored in subsequent ramp-fitting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "jump_map.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "group_indexes = np.arange(jump_map.shape[1]).astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pick one pixel with a flagged jump, and find the group(s) with the jump flags\n",
    "j_index = 10112\n",
    "jumpy = jump_map_indexes[0][j_index]\n",
    "jumpx = jump_map_indexes[1][j_index]\n",
    "jump_grp = jump_map[0, :, jumpy, jumpx]\n",
    "print('Jump located in group(s) {} of pixel ({}, {})'.format(group_indexes[jump_grp], jumpx, jumpy))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the signal up the ramp for this pixel\n",
    "plot_jump(jump.data[0, :, jumpy, jumpx], jump_grp, xpixel=jumpx, ypixel=jumpy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot a grid of some examples of flagged jumps. The red marks signify flagged jumps. These signal values will be ignored in subsequent ramp-fitting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "indexes_to_plot = [200, 401, 600, 30010, 31000, 1202, 1400, 21600, 10112]\n",
    "jump_data = np.zeros((jump.shape[1], len(indexes_to_plot)))\n",
    "jump_grps = np.zeros((jump.shape[1], len(indexes_to_plot))).astype(bool)\n",
    "jump_locs = []\n",
    "for counter, idx in enumerate(indexes_to_plot):\n",
    "    #integ, grp, y, x = jump_flags[idx]\n",
    "    y = jump_map_indexes[0][idx]\n",
    "    x = jump_map_indexes[1][idx]\n",
    "    grp = jump_map[0, :, y, x]\n",
    "\n",
    "    jump_data[:, counter] = jump.data[0, :, y, x]\n",
    "    jump_grps[:, counter] = grp\n",
    "    jump_locs.append((x, y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plot_jumps(jump_data, jump_grps, jump_locs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Top of Notebook](#top)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='ramp_fitting'></a>\n",
    "## The `Ramp Fitting` step"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Summary\n",
    "\n",
    "This step performs line-fitting to the corrected data, and produces a slope image for each integration. For a given pixel, any groups within an integration that contain a jump flag or that are flagged as saturated are ignored.\n",
    "\n",
    "For the purposes of this notebook, we will use the `save_opt` parameter to tell the ramp-fitting step to save an optional output file that contains some of the details on the ramp fits. This information will be used for the plots after the step is run. By default, `save_opt` is False and the optional outputs are not saved.\n",
    "\n",
    "\n",
    "#### Documentation\n",
    "\n",
    "[Full description](https://jwst-pipeline.readthedocs.io/en/stable/jwst/ramp_fitting/description.html) of the step.\n",
    "\n",
    "#### Arguments\n",
    "\n",
    "The jump step has [several optional arguments](https://jwst-pipeline.readthedocs.io/en/stable/jwst/ramp_fitting/arguments.html), including the ability to save optional outputs into a second file.\n",
    "\n",
    "#### Reference files used\n",
    "\n",
    "This step uses the [`READNOISE`](https://jwst-pipeline.readthedocs.io/en/stable/jwst/references_general/readnoise_reffile.html) and [`GAIN`](https://jwst-pipeline.readthedocs.io/en/stable/jwst/references_general/gain_reffile.html) reference files. These files contain maps of the readnoise and gain values across the detector.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Run the step"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "List the available parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(RampFitStep.spec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using the run() method\n",
    "ramp_fit_step = RampFitStep()\n",
    "ramp_fit_step.output_dir = output_dir\n",
    "ramp_fit_step.save_results = True\n",
    "\n",
    "# Let's save the optional outputs, in order\n",
    "# to help with visualization later\n",
    "ramp_fit_step.save_opt = True\n",
    "\n",
    "# Call using the dark instance from the previously-run\n",
    "# jump step\n",
    "ramp_fit = ramp_fit_step.run(jump)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An important detail here is that there are two output files: \n",
    "\n",
    "1. The file ending with `*_0_rampfitstep.fits` contains the mean slope image across all integrations in the exposure.\n",
    "\n",
    "2. The file ending with `*_1_rampfitstep.fits` contains a separate slope image for each integration. In this case our exposure contains only a single integration, so the data in the two files are identical."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rampfit_output_file = os.path.join(output_dir, '{}_0_rampfitstep.fits'.format(input_file_base))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are working with the output datamodels in this notebook. Unlike the preceeding steps, where the output was a single datamodel, the ramp_fitting step outputs a tuple of 2 datamodel instances. The first element in the tuple is the datamodel instance containing the mean rate image, while the second element is the datamodel instance containing a separate slope image for each integration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(ramp_fit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(ramp_fit[0]), type(ramp_fit[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Show the shape of the two output data models. The output with the mean slope image is only 2-dimensional, while the output with one slope image for each integration is 3-dimensional, even in this case where we have only one integration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ramp_fit[0].shape, ramp_fit[1].shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we see that in our case with a single integration, the two datamodel instances contain identical data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ramp_fit[0].data[500, 500], ramp_fit[1].data[0, 500, 500]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And note that the data in these files are the slopes from the ramp fitting, so the units have changed from DN to DN/sec."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show the units of the ramp-fit data\n",
    "ramp_fit[0].meta.bunit_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's load the data from the third output file, which contains the optional outputs. Our goal is to grab the intercept values from the line-fitting and use them to re-create the plots from the jump step and overplot the best-fit signals on top of the signals that went into the step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "optional_file = os.path.join(output_dir, '{}fitopt.fits'.format(input_file_base))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hdulist = fits.open(optional_file)\n",
    "hdulist.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The intercepts from the line-fitting are stored in the `YINT` extension."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hdulist['YINT'].data[0, :, 200, 200]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For our plots below, we need the intercept values from the line-fitting. Let's ignore all but the first plane of the extension, since those are zeros."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "intercepts = hdulist['YINT'].data[0, 0, :, :]\n",
    "hdulist.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the pixels to be plotted, create linear ramps using the output slopes and intercepts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the exposure time associated with each group\n",
    "num_groups = ramp_fit[0].meta.exposure.ngroups\n",
    "group_time = ramp_fit[0].meta.exposure.group_time\n",
    "group_times = np.arange(num_groups) * group_time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reconstruct linear ramps from the slope and intercept values\n",
    "lin_ramps = np.zeros((jump.shape[1], len(indexes_to_plot)))\n",
    "for counter, idx in enumerate(indexes_to_plot):\n",
    "    y = jump_map_indexes[0][idx]\n",
    "    x = jump_map_indexes[1][idx]\n",
    "    grp = jump_map[0, :, y, x]\n",
    "\n",
    "    rate = ramp_fit[0].data[y, x]\n",
    "    intercept = intercepts[y, x]\n",
    "    lin_ramps[:, counter] = intercept + (rate * group_times)    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lin_data = intercepts[jumpy, jumpx] + (ramp_fit[0].data[jumpy, jumpx] * group_times)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_jump(jump.data[0, :, jumpy, jumpx], jump_grp, xpixel=jumpx,\n",
    "          ypixel=jumpy, slope=lin_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_jumps(jump_data, jump_grps, jump_locs, slopes=lin_ramps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at the slope image. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_image(ramp_fit[0].data, 0, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What's going on with the dark pixels scattered across the left side of the detector? It turns out they have signal values of exactly zero in the slope image. Let's have a closer look, and ignore the reference pixels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scipix = ramp_fit[0].data[4:2044, 4:2044]\n",
    "zero_pix = np.where(scipix == 0.0)\n",
    "print('{} science pixels have a slope that is exactly zero.'.format(len(zero_pix[0])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pick one of these pixels, and examine the pixel's signals up the ramp, as well as its data quality flags up the ramp."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = 45\n",
    "\n",
    "# Add 4 to the coordinates because we stripped off the 4 columns/rows of \n",
    "# refpix in the np.where statement above\n",
    "y = 4 + zero_pix[0][idx]\n",
    "x = 4 + zero_pix[1][idx]\n",
    "print('({}, {}) has a slope of 0.'.format(x, y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print the slope (should be zero), signal values up the ramp, and DQ flags up the ramp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scipix[y-4, x-4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "linearity.data[0, :, y, x]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "linearity.groupdq[0, :, y, x]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What does a DQ flag value of 2 mean?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dqflags.pixel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So this pixel (and the others with values of 0 in the slope image) have been flagged as saturated in all groups. In this case, the pipeline cannot calculate a slope value, and assigns a value of zero. We'll see in the Stage 3 notebook how these pixels will be ignored when combining multiple exposures to create a final mosaic image."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Top of Notebook](#top)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bonus Topic: Logging"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We've seen in the calls to the pipeline and steps above that the collection of logging entries output to the screen can be quite long. If you wish to capture this information in a file rather than on the screen, then you can create a simple logging configuration file that resides in the directory where you are calling the pipeline. By setting the `logcfg` keyword to the name of this file in your pipeline calls, the logging entries will be routed to the file specified in the configuration file.\n",
    "\n",
    "See the [log file set up instructions](https://jwst-pipeline.readthedocs.io/en/latest/jwst/introduction.html#logging-configuration) for details."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='exercise_solutions'></a>\n",
    "## Exercise Solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='exercise1_solution'></a>\n",
    "### Exercise 1: Run MIRI data through the pipeline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the supplied MIRI exposure (which can be downloaded in the [Download Data](#download_data) section), run the calwebb_detector1 pipeline, examnine which steps are run, and note the different order compared to the pipeline run above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Using the run() method:\n",
    "# Instantiate the pipeline\n",
    "miri1 = calwebb_detector1.Detector1Pipeline()\n",
    "\n",
    "# Save the final output of the pipeline\n",
    "miri1.output_dir = output_dir\n",
    "miri1.save_results = True\n",
    "\n",
    "# Save the output from the dark current subtraction step,\n",
    "# and let's use a more stringent limit for cosmic ray\n",
    "# detection\n",
    "miri1.dark_current.save_results = True\n",
    "miri1.jump.rejection_threshold = 3\n",
    "\n",
    "# Call the run() method\n",
    "miri_output = miri1.run(miri_uncal_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using the call() method:\n",
    "# In order to match the parameters used in the the call to run() above,\n",
    "# you first must edit the parameter reference file, which is:\n",
    "miri_det1_param_reffile"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add the following lines to the bottom of the file:\n",
    "\n",
    "```\n",
    "- class: jwst.dark_current.dark_current_step.DarkCurrentStep\n",
    "  name: dark_current\n",
    "  parameters: {save_results: True}\n",
    "- class: jwst.jump.jump_step.JumpStep\n",
    "  name: jump\n",
    "  parameters: {rejection_threshold: 3}\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now use the call() method and supply the parameter reference file\n",
    "miri_output = calwebb_detector1.Detector1Pipeline.call(miri_uncal_file,\n",
    "                                                       output_dir=output_dir,\n",
    "                                                       save_results=True,\n",
    "                                                       config_file=miri_det1_param_reffile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Look at one of the output rate files\n",
    "# HINT: Open one of the rate files and use show_image with\n",
    "# min and max signal rate limits of 0 and 100.\n",
    "miri_rate_file = miri_uncal_file.replace('uncal.fits', 'rate.fits')\n",
    "miri_rate = datamodels.open(miri_rate_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_image(miri_rate.data, 0, 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Back to [Exercise 1](#exercise1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='exercise2_solution'></a>\n",
    "### Exercise 2: Re-run the `reference pixel subtraction` step"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try re-running the [reference pixel subtraction](#refpix) step and changing the `side_smoothing_length` and/or `side_gain` parameters. Examine the output to see how well the 1/f noise is removed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check the names of the available parameters\n",
    "RefPixStep.spec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using the run() method.\n",
    "\n",
    "# Instantiate the step\n",
    "refpix_step = RefPixStep()\n",
    "\n",
    "# Set the side_smoothing_length to a non-default value\n",
    "# (default is 11)\n",
    "refpix_step.side_smoothing_length = 50\n",
    "\n",
    "# Call using the superbias instance from the previously-run\n",
    "# saturation step\n",
    "refpix_ex2 = refpix_step.run(superbias)\n",
    "# OR run using the output file from the superbias step\n",
    "# refpix_ex2 = refpix_step.run(superbias_output_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using the call() method\n",
    "\n",
    "# First open 'refpix_paramfile.asdf' and set the \n",
    "# side_smoothing_length entry to a new value.\n",
    "refpix_ex2 = RefPixStep.call(superbias, output_dir=output_dir,\n",
    "                                save_results=True, \n",
    "                                config_file=refpix_param_reffile)\n",
    "\n",
    "# OR run using the output file from the superbias step\n",
    "#refpix_ex2 = RefPixStep.call(superbias_output_File, output_dir=output_dir,\n",
    "#                                save_results=True, \n",
    "#                                config_file=refpix_param_reffile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Back to [Exercise 2](#exercise2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That is the end of the Stage 1 pipeline. We have now transformed a single exposure from a raw, multiaccum ramp into a signal rate image after applying detector-level corrections. From here, the data move on to the Stage 2 pipeline. For imaging data such as these, that means the *calwebb_image2* pipeline. This will be shown in the Stage 2 notebook. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Top of Notebook](#top)"
   ]
  }
 ],
 "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}