{ "cells": [ { "cell_type": "markdown", "id": "3e21d550", "metadata": {}, "source": [ "\"stsci_logo\" \n", "\n", "# TSO JWebbinar Notebook 2: Exploring TSO Products\n", "-----\n", "\n", "**Author**: Néstor Espinoza, AURA Assistant Astronomer, NIRISS branch\n", "
\n", "**Pipeline Version**: 1.3.3 on March 2, 2022\n", "
\n", "**Last Updated**: version 1.5.2 on June 16, 2022" ] }, { "cell_type": "markdown", "id": "acc24f5c", "metadata": {}, "source": [ "## Table of contents\n", "1. [Introduction](#intro)
\n", " 1.1 [Purpose of this Notebook](#purpose)
\n", " 1.2 [Input Simulations](#inputs)
\n", " 1.3 [Caveats for Simulated Data](#nirisscaveats)
\n", "2. [Setup](#setup)
\n", "3. [A TSO tour through the `Detector1` stage](#detector1)
\n", " 3.1 [Checking data-quality flags](#dqflags)
\n", " 3.2 [Reference pixel corrections](#refpix)
\n", " 3.3 [Linearity corrections](#linearity)
\n", " 3.4 [1/f noise](#oneoverf)
\n", " 3.5 [The \"jump\" detection](#jump)
\n", " 3.6 [The ramp-fitting step](#rampfit)
\n", "4. [A (quick) TSO look at the `Spec2` stage](#spec2)
\n", "5. [General TSO product exploration](#tsoexp)
\n", " 5.1 [Spectral tracing](#spectracing)
\n", " 5.2 [Coming back to 1/f noise](#one-over-f2)
\n", " 5.3 [Spectral extraction](#specextraction)
\n", " 5.4 [Quick lightcurve exploration](#lcexploration)
\n", " 5.5 [Quick lightcurve sensitivity exploration](#sensexploration)
\n", " 5.6 [White-light lightcurve residual analyses](#lcresiduals)
\n", "6. [Final words](#final-words)
" ] }, { "cell_type": "markdown", "id": "0c0a705b", "metadata": {}, "source": [ "---\n", "\n", "1.-Introduction \n", "------------------" ] }, { "cell_type": "markdown", "id": "1ffeef35", "metadata": {}, "source": [ "### 1.1-Purpose of this Notebook ###\n", "\n", "In this notebook we provide material to get started on the exploration of Time-Series Observations (TSO) products; in particular, we focus on outputs from `detector1` stage of the pipeline, and go through our own spectral tracing and extraction procedures, in order to arrive to general lightcurve analyses JWST TSO users might want to perform on their datasets." ] }, { "cell_type": "markdown", "id": "cfc00ab4", "metadata": {}, "source": [ "### 1.2-Input Simulations ###\n", "\n", "The input data for this notebook are simulations made by the NIRISS/SOSS team of a transit of WASP-43b. These simulations were made using the [`mirage`](https://mirage-data-simulator.readthedocs.io/en/latest/) simulator. This full dataset can be found [here](https://stsci.box.com/s/fbolnundf2q01yq3fkqpbvswhqch7c91), and includes both the `uncal` products (i.e., the \"raw\" data) and the data calibrated after the `assignwcs` step of the [JWST pipeline](https://jwst-pipeline.readthedocs.io/en/latest/jwst/pipeline/calwebb_spec2.html). Note how there are two datasets for each of those set of products, as the JWST products are divided into \"segments\" dividing the entire exposure into smaller chunks for simple data size purposes.\n", "\n", "We will be using products at two different stages of the pipeline, namely:\n", "\n", "1. **The Uncalibrated Products**. These correspond to data that have not yet been passed through any of the stages of the pipeline. These have `*uncal.fits` extensions.\n", "3. **The partially Stage 2 Calibrated Products**. These correspond to data that have only been calibrated up to the `assign_wcs` step, with which one can obtain the wavelength map of the spectra. These have `*_1_assignwcsstep.fits` extensions.\n", "\n", "These products are downloaded directly from this notebook below. Importantly, these simulations were generated with the following reference files (we will define what these are and why they are important to define below!):\n", "\n", "- Reference file for the `SuperBias` step: `jwst_niriss_superbias_0017.fits`.\n", "\n", "- Reference file for the `Linearity` step: `jwst_niriss_linearity_0013.fits`." ] }, { "cell_type": "markdown", "id": "e82389a0", "metadata": {}, "source": [ "### 1.3-Caveats for Simulated Data ###\n", "\n", "The [`mirage`](https://mirage-data-simulator.readthedocs.io/en/latest/) simulator is not perfect. Be aware that how close these simulations are to real data is not going to be perfectly known until we obtain on-sky data with JWST. As such, a big disclaimer is that while `mirage` simulations represent a large fraction of the systematics that we expect to see, it is in no way a guarantee that this will be the full picture of what we will _actually_ see on real JWST data.\n", "\n", "Also note that time-series systematics have not been added to those simulations. This has been important for previous observatories to consider, and appear due to, e.g., thermal breathing of the instruments, focus changes, etc. " ] }, { "cell_type": "markdown", "id": "62f3de65", "metadata": {}, "source": [ "2.-Setup \n", "------------------\n", "\n", "In this section we set a number of necessary things in order for the pipeline to run successfully.\n", "\n", "First we'll set the CRDS context; this dictates the versions of various pipeline reference files to use. \n", "\n", "Next we'll import the various python packages that we're actually going to use in this notebook, including both generic utility functions and the actual pipeline modules themselves. \n", "\n", "Finally, we also define two directories: the `working_folder` (where we will be saving the outputs in this notebook) and our `data_folder`, where we can have some pre-loaded data products. If you want to run this notebook from scratch, those two folders should be the same --- this is how it is right now by default. If you have some pre-loaded products (e.g., you were part of the JWebbinar and you have access to pre-loaded outputs), then you can define the `data_folder` to point to where those files are:" ] }, { "cell_type": "code", "execution_count": null, "id": "c084f9b6", "metadata": {}, "outputs": [], "source": [ "# Need to set these enviromental variables for this notebook to work properly:\n", "%set_env CRDS_PATH /path/to/my/folder/\n", "%set_env CRDS_SERVER_URL https://jwst-crds-pub.stsci.edu\n", " \n", "# Import all useful libraries:\n", "import os\n", "import pickle\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from numpy.polynomial import chebyshev\n", "from scipy.ndimage import median_filter\n", "from scipy.ndimage import gaussian_filter1d\n", "from mpl_toolkits.axes_grid1.inset_locator import mark_inset\n", "from mpl_toolkits.axes_grid1 import make_axes_locatable\n", "from scipy.interpolate import interp1d, splrep,splev\n", "\n", "from jwst import datamodels\n", "from astropy.io import fits\n", "from astropy.timeseries import LombScargle\n", "from astropy.utils.data import download_file\n", "\n", "# Import the JWST pipeline stages (which contain the steps):\n", "from jwst.pipeline import calwebb_detector1\n", "from jwst.pipeline import calwebb_spec2\n", "\n", "# Define the folder on which we will be working on (here we'll store all our outputs):\n", "working_folder = os.getcwd()+'/'\n", "data_folder = os.getcwd()+'/'" ] }, { "cell_type": "markdown", "id": "d49194a6", "metadata": {}, "source": [ "3.-A TSO tour through the `Detector1` stage \n", "------------------\n", "\n", "Here, we will start reducing data and checking the results of this data reduction for TSO's on each step. This will allow us to pick our brains as to how to explore these TSO products, but also, to be careful to be mindful of what the pipeline is actually doing --- and how impactful/important this is for TSOs. \n", "\n", "Here, we will work with only one segment of the data --- on a real data analysis session, we would need to run the same steps on both segments to get the reductions for the entire exposure. We download/load the uncalibrated second segment of this exposure below:" ] }, { "cell_type": "code", "execution_count": null, "id": "e5a0097f", "metadata": {}, "outputs": [], "source": [ "# Download uncal file:\n", "file_link = 'https://stsci.box.com/shared/static/dy0t8hza8tzpoqcdlny54e9y6pbsjctc.fits'\n", "filename = 'WASP43_NIS_SOSS-seg002_CLEAR_uncal.fits'\n", "\n", "# Look if the file has been downloaded. If not, do it:\n", "if os.path.exists(data_folder + filename):\n", " \n", " print(filename + ' already exists. Not downloading.')\n", "\n", "else:\n", " \n", " # Download file:\n", " print('Downloading {}...'.format(filename))\n", " download_filename = download_file(file_link, cache=False)\n", " \n", " # Rename file:\n", " os.rename(download_filename, data_folder + '/' + filename)" ] }, { "cell_type": "markdown", "id": "1530410f", "metadata": {}, "source": [ "Done! Now let's load the data inside a JWST `datamodel`:" ] }, { "cell_type": "code", "execution_count": null, "id": "57b4af42", "metadata": {}, "outputs": [], "source": [ "uncal_data = datamodels.RampModel(data_folder + '/' + filename)" ] }, { "cell_type": "code", "execution_count": null, "id": "f190a25e", "metadata": {}, "outputs": [], "source": [ "# Cover a bug that exists given the simulated data was created a while ago. This will not be \n", "# needed to be done on real data, but it's a nice excercise showing how easy it is to deal \n", "# with JWST datamodels:\n", "uncal_data.meta.dither.dither_points = int(uncal_data.meta.dither.dither_points)" ] }, { "cell_type": "markdown", "id": "4be76375", "metadata": {}, "source": [ "## 3.1 Checking data-quality flags\n", "\n", "An important component of any TSO analysis is to flag bad pixels, pixels identified as cosmic rays and/or identify saturated pixels. The former pixels are actually ingested directly by instrument experts on a mask that one loads on the very first step of the JWST pipeline for NIRISS/SOSS data: the `dq_init_step`. Let's run it:" ] }, { "cell_type": "code", "execution_count": null, "id": "f0a184ce", "metadata": {}, "outputs": [], "source": [ "# Let's run the DQ init step --- this will assign DQ flags to our fits uncal file. \n", "# First, create output directory:\n", "\n", "if not os.path.exists(working_folder + '/' + 'pipeline_outputs_directory'):\n", " os.mkdir(working_folder + 'pipeline_outputs_directory')\n", " \n", "pipeline_outputs_folder = working_folder + 'pipeline_outputs_directory'\n", "\n", "# Run step:\n", "dq_results = calwebb_detector1.dq_init_step.DQInitStep.call(uncal_data, \n", " output_dir=pipeline_outputs_folder, \n", " save_results=True)" ] }, { "cell_type": "markdown", "id": "ebfcb7cb", "metadata": {}, "source": [ "The full list of data quality flags is [here](https://jwst-pipeline.readthedocs.io/en/latest/jwst/references_general/references_general.html?highlight=data%20quality%20flags#data-quality-flags). For example, pixels with data quality flags of `0` are \"normal\" pixels and pixels with data quality flags of `2147483648` are reference pixels (more on those in a bit). For all the rest, one has to have some special care. They might be used, but one might want to be careful on how to move forward with them on a case by case basis. \n", "\n", "### 3.1.1 Checking \"special\" pixels\n", "\n", "**Let's assume we want to be extremely conservative, and plot all pixels which are _not_ normal (i.e., having a value other than `0`) or reference pixels (i.e., having a value different than `2147483648`)**. This can be done by extracting the pixel data-quality flags, and using their values as follows: " ] }, { "cell_type": "code", "execution_count": null, "id": "56e8c6e1", "metadata": {}, "outputs": [], "source": [ "# To get the bad pixels, we need to use the \"pixeldq\" which sets the general data quality for each pixel:\n", "dq_pixels = dq_results.pixeldq\n", "\n", "# Identify pixels which are NOT OK and are NOT reference pixels --- we call those \"special\" pixels:\n", "idx_special_pixels = np.where((dq_pixels != 0)&(dq_pixels != 2147483648))" ] }, { "cell_type": "markdown", "id": "960fe4b3", "metadata": {}, "source": [ "Let's plot where these are located in the frames:" ] }, { "cell_type": "code", "execution_count": null, "id": "44706803", "metadata": {}, "outputs": [], "source": [ "# Let's plot these \"special\" pixels now:\n", "plt.figure(figsize=(20,5))\n", "\n", "# Extract a zero-array with the same shape as the original frames; use first group, first integration as reference:\n", "special_pixels = np.zeros(dq_results.data[0, 0, :, :].shape)\n", "\n", "# Identify \"special\" pixels as ones in this array (initially) zero-array:\n", "special_pixels[idx_special_pixels] = 1.\n", "\n", "# Plot them:\n", "im = plt.imshow(special_pixels, interpolation = 'None')\n", "im.set_clim(0,1)\n", "plt.title('Special pixels in our NIRISS/SOSS SUBSTRIP256 frame:')" ] }, { "cell_type": "markdown", "id": "9e23c07a", "metadata": {}, "source": [ "As can be seen, most of their location are random. There seems to be a set of pixels in the corners that define the boundary with the reference pixels, however, which one should be wary of.\n", "\n", "### 3.1.1 Checking saturated pixels\n", "\n", "One very important detail in JWST data analysis involves checking which pixels are \"saturated\" or not. Saturation in the JWST context is an [instrument-by-instrument defined upper signal level](https://jwst-docs.stsci.edu/methods-and-roadmaps/jwst-time-series-observations/tso-saturation) --- but is nonetheless good to know that one can control what is considered saturated or not, and also to check which pixels are identified as saturateed by the pipeline. This is done through the `saturation` step --- the next step of the JWST pipeline for Detector 1. Let's run it:" ] }, { "cell_type": "code", "execution_count": null, "id": "00464d78", "metadata": {}, "outputs": [], "source": [ "# Run step:\n", "saturation_results = calwebb_detector1.saturation_step.SaturationStep.call(dq_results , \n", " output_dir=pipeline_outputs_folder, \n", " save_results=True)" ] }, { "cell_type": "markdown", "id": "a75bfef9", "metadata": {}, "source": [ "The saturation step works by comparing the observed ADU values with the saturation signal-levels defined for each pixel in a reference file. As can be seen above, that reference file is indicated by the line `stpipe.SaturationStep - INFO - Using SATURATION reference file [yourfile]`. In the case of our run at the time of writing, this was the `jwst_niriss_saturation_0011.fits ` file --- but this might change during commisioning or at later stages of pipeline/instrument development and analysis.\n", "\n", "Let's have a look at what pixels are indicated **only** as saturated in the first and last groups. This is obtained through the `groupdq` product (which is different from the `pixeldq` used above) --- which are data-quality flags _per group_. Values of `2` mean the pixel is saturated in a given group:" ] }, { "cell_type": "code", "execution_count": null, "id": "17085a32", "metadata": {}, "outputs": [], "source": [ "# Extract group data-quality:\n", "dq_groups = saturation_results.groupdq\n", "\n", "# Identify saturated pixels in the first and last groups:\n", "idx_sat_pixels_first_group = np.where(dq_groups[0, 0, :, :] == 2)\n", "idx_sat_pixels_last_group = np.where(dq_groups[0, -1, :, :] == 2)" ] }, { "cell_type": "markdown", "id": "5b0a4728", "metadata": {}, "source": [ "Let's check their location:" ] }, { "cell_type": "code", "execution_count": null, "id": "3dd84b2b", "metadata": {}, "outputs": [], "source": [ "# Extract a zero-array with the same shape as the original frames; use first group, first integration as reference:\n", "saturated_pixels_firstgroup = np.zeros(saturation_results.data[0, 0, :, :].shape)\n", "saturated_pixels_lastgroup = np.zeros(saturation_results.data[0, 0, :, :].shape)\n", "\n", "# Identify \"saturated\" pixels as ones in this array (initially) zero-array:\n", "saturated_pixels_firstgroup[idx_sat_pixels_first_group] = 1.\n", "\n", "# Plot them:\n", "plt.figure(figsize=(20,5))\n", "im = plt.imshow(saturated_pixels_firstgroup, interpolation = 'None')\n", "im.set_clim(0,1)\n", "plt.title('Saturated pixels on FIRST group, first integration:')\n", "\n", "# Repeat for last group:\n", "saturated_pixels_lastgroup[idx_sat_pixels_last_group] = 1.\n", "\n", "# Plot them:\n", "plt.figure(figsize=(20,5))\n", "im = plt.imshow(saturated_pixels_lastgroup, interpolation = 'None')\n", "im.set_clim(0,1)\n", "plt.title('Saturated pixels on LAST group, first integration:')" ] }, { "cell_type": "markdown", "id": "ed6d8613", "metadata": {}, "source": [ "As can be seen, there are not many pixels saturated in the first group --- and pixels saturated in the first group remain saturated in the last group. This is the case for the pixel close to row 250 and column 500.\n", "\n", "Now --- why were the pixels identified as saturated? This comes from the `saturation` reference file. You can check which one was used directly from the outputs of the pipeline as follows:\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "cd3b1436", "metadata": {}, "outputs": [], "source": [ "saturation_results.meta.ref_file.saturation.name" ] }, { "cell_type": "markdown", "id": "4e277e99", "metadata": {}, "source": [ "What does this file looks like? Let's check it out:" ] }, { "cell_type": "code", "execution_count": null, "id": "42ae83fc", "metadata": {}, "outputs": [], "source": [ "# Base directory where reference files are stored (this was defined in the Setup section above:)\n", "# Add /references/jwst/niriss/\n", "base_ref_files = '/path/to/my/folder/references/jwst/niriss/'\n", "\n", "# Read it in:\n", "saturation_ref_file = fits.open(base_ref_files+saturation_results.meta.ref_file.saturation.name[7:])" ] }, { "cell_type": "code", "execution_count": null, "id": "91279ea2", "metadata": {}, "outputs": [], "source": [ "# Let's print values of the reference file:\n", "print('File values:')\n", "print(saturation_ref_file[1].data[:,:])\n", "\n", "print('File shape:')\n", "print(saturation_ref_file[1].data.shape)" ] }, { "cell_type": "markdown", "id": "ba32975e", "metadata": {}, "source": [ "As you see above, the reference file is a file that, for each pixel, has a value that marks when it's \"saturated\". This value seems to be around 70,000 counts. Note also that the shape is (2048, 2048) --- the saturation values are defined by the entire full frame of the NIRISS detectors --- this file is simply \"cut\" to match the dimensions of the subarray (which is a subset of the pixels of the full frame) to define the saturation values for the `SUBSTRIP256` subarray. Let's make this cut directly:" ] }, { "cell_type": "code", "execution_count": null, "id": "0d1d409c", "metadata": {}, "outputs": [], "source": [ "fig, ax1 = plt.subplots(figsize=(20, 5))\n", "\n", "im1 = ax1.imshow(saturation_ref_file[1].data[-256:,:], cmap='Reds', vmin = 50000, vmax = 70000, interpolation = None)\n", "\n", "# Colorbar:\n", "divider = make_axes_locatable(ax1)\n", "cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n", "\n", "cbar = fig.colorbar(im1, shrink = 0.08, cax=cax)\n", "cbar.ax.tick_params(labelsize=13)\n", "cbar.ax.get_yaxis().labelpad = 20\n", "cbar.ax.set_ylabel('Saturation counts', rotation=270, fontsize = 15)" ] }, { "cell_type": "markdown", "id": "ddf46ae5", "metadata": {}, "source": [ "Lots of interesting structure here! Most of it follows the bias (as we'll see below), but there are also \"bubbles\" of higher saturation values in the image. Reference pixels (in the left, right and bottom corners) also saturate at larger flux values.\n", "\n", "Saturation is important because later, when running the ramp-fitting step, saturated groups are ommited from the calculations. You can, of course, manually modify the reference file and define lower count levels for saturation if you want. Alternatively, you can also modify the `groupdq` values directly of the output objects the pipeline gives --- when you ingest them to later steps of the pipeline, they will be propagated accordingly. \n", "\n", "For instance, we could mark all the groups above the third of pixel `(123,123)` as saturated in the first integration as follows:" ] }, { "cell_type": "code", "execution_count": null, "id": "566c15bd", "metadata": {}, "outputs": [], "source": [ "saturation_results.groupdq[0, 2:, 123, 123] = 2" ] }, { "cell_type": "markdown", "id": "0c65ab2c", "metadata": {}, "source": [ "**3.1.1 Exercise for the reader: above, we identified a handful of pixels that were marked as saturated by the pipeline. Go back to those pixels, and check that their actual values at the groups where they are identified as saturated indeed are larger than what the saturation reference file defines as \"saturated\".**\n", "\n", "**3.1.2 Exercise for the reader: we only checked pixel values that are saturated above. However, what if they are saturated _and_ they have other data quality flags as well? Read the documentation and find out why a call such as `np.where((saturation_results.groupdq[0, :, :, :] & 4) > 0)` would be a more generic call to identify _all_ saturated pixels.**" ] }, { "cell_type": "markdown", "id": "bf0355b1", "metadata": {}, "source": [ "## 3.2 Reference pixel corrections\n", "\n", "After the saturation step, there are two steps to go to get to the step we want to get at --- the `refpix` step. \n", "We first perform the `superbias` step which removes the detector pedestal from our groups. We are careful to consider that, as it was said in the introduction, the superbias reference file is a special one: the `jwst_niriss_superbias_0017.fits` reference file. We have uploaded this reference file to Box, but this can be obtained directly from CRDS:" ] }, { "cell_type": "code", "execution_count": null, "id": "fc23ba69", "metadata": {}, "outputs": [], "source": [ "# Download superbias reference file:\n", "ref_filename = 'jwst_niriss_superbias_0017.fits'\n", "download_link = 'https://stsci.box.com/shared/static/1n3etewqfhgbdhd1vpwour3gg3shiisa.fits'\n", "\n", "# Look if the file has been downloaded. If not, do it:\n", "if os.path.exists(data_folder + ref_filename):\n", " \n", " print(ref_filename + ' already exists. Not downloading.')\n", "\n", "else:\n", " \n", " print('Downloading {}...'.format(ref_filename)+' from '+ download_link)\n", " download_filename = download_file(download_link,\n", " cache=False)\n", " # Rename file:\n", " os.rename(download_filename, data_folder + ref_filename)" ] }, { "cell_type": "code", "execution_count": null, "id": "47e0fe3e", "metadata": {}, "outputs": [], "source": [ "# Run superbias step:\n", "superbias_results = calwebb_detector1.superbias_step.SuperBiasStep.call(saturation_results , \n", " output_dir=pipeline_outputs_folder, \n", " override_superbias = data_folder + ref_filename,\n", " save_results=True)" ] }, { "cell_type": "markdown", "id": "73e8915e", "metadata": {}, "source": [ "Let's check how different the superbias and the saturation products look. To this end, let's plot the last group of the first integration:" ] }, { "cell_type": "code", "execution_count": null, "id": "079d7caa", "metadata": {}, "outputs": [], "source": [ "# Plot them:\n", "plt.figure(figsize=(20,5))\n", "im = plt.imshow(saturation_results.data[0,-1,:,:] / np.nanmedian(saturation_results.data[0,-1,:,:]), \\\n", " interpolation = 'None')\n", "im.set_clim(-3,2)\n", "plt.title('Before the Superbias step:')\n", "\n", "# Repeat for last group:\n", "saturated_pixels_lastgroup[idx_sat_pixels_last_group] = 1.\n", "\n", "# Plot them:\n", "plt.figure(figsize=(20,5))\n", "im = plt.imshow(superbias_results.data[0,-1,:,:] / np.nanmedian(superbias_results.data[0,-1,:,:]), \\\n", " interpolation = 'None')\n", "im.set_clim(-3,2)\n", "plt.title('After the Superbias step:')" ] }, { "cell_type": "markdown", "id": "340e93a6", "metadata": {}, "source": [ "Wow! That's a huge change. Overall, this looks much better. Let's plot the profiles the column whose index is 1500 to have a closer look:" ] }, { "cell_type": "code", "execution_count": null, "id": "85918a2a", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(nrows = 1, ncols = 2, figsize=(15,5))\n", "\n", "ax[0].plot(saturation_results.data[0,-1,:,1500], label = 'Before the Superbias step')\n", "ax[0].plot(superbias_results.data[0,-1,:,1500], label = 'After the Superbias step')\n", "ax[0].set_xlabel('Row pixel index', fontsize = 14)\n", "ax[0].set_ylabel('Counts', fontsize = 14)\n", "ax[0].set_title('Comparison before/after Superbias step', fontsize = 14)\n", "ax[0].legend()\n", "\n", "ax[1].set_title('Same, but median-substracted counts', fontsize = 14)\n", "ax[1].plot(saturation_results.data[0,-1,:,1500] - np.nanmedian(saturation_results.data[0,-1,:,1500]))\n", "ax[1].plot(superbias_results.data[0,-1,:,1500] - np.nanmedian(superbias_results.data[0,-1,:,1500]))\n", "ax[1].set_xlabel('Row pixel index', fontsize = 14)\n", "ax[1].set_ylabel('Counts - Median Counts', fontsize = 14)" ] }, { "cell_type": "markdown", "id": "4d5ba157", "metadata": {}, "source": [ "As can be seen, a ton of structure has been removed. Also, all the pixels seem to be at the same background level. This is a good sign that the Superbias correction has worked correctly. Note the background level is _below_ zero, but this is OK --- group-to-group changes in the overall flux level _are_ expected, and will be removed in the next step --- the `ref_pix` step. Let's apply it and see what comes out of this:" ] }, { "cell_type": "code", "execution_count": null, "id": "9b133a2b", "metadata": {}, "outputs": [], "source": [ "# Run superbias step:\n", "refpix_results = calwebb_detector1.refpix_step.RefPixStep.call(superbias_results, \n", " output_dir=pipeline_outputs_folder, \n", " save_results=True)" ] }, { "cell_type": "markdown", "id": "67c3e82f", "metadata": {}, "source": [ "Let's do the same comparison we did above after the Superbias step:" ] }, { "cell_type": "code", "execution_count": null, "id": "ad093d19", "metadata": {}, "outputs": [], "source": [ "# Plot them:\n", "plt.figure(figsize=(20,5))\n", "im = plt.imshow(superbias_results.data[0,-1,:,:] / np.nanmedian(superbias_results.data[0,-1,:,:]), \\\n", " interpolation = 'None')\n", "im.set_clim(-3,2)\n", "plt.title('Before the RefPix step:')\n", "\n", "# Plot them:\n", "plt.figure(figsize=(20,5))\n", "im = plt.imshow(refpix_results.data[0,-1,:,:] / np.nanmedian(refpix_results.data[0,-1,:,:]), \\\n", " interpolation = 'None')\n", "im.set_clim(-3,2)\n", "plt.title('After the RefPix step:')" ] }, { "cell_type": "code", "execution_count": null, "id": "7757ed99", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(nrows = 1, ncols = 2, figsize=(15,5))\n", "\n", "ax[0].plot(superbias_results.data[0,-1,:,1500], label = 'Before the RefPix step')\n", "ax[0].plot(refpix_results.data[0,-1,:,1500], label = 'After the RefPix step')\n", "ax[0].set_xlabel('Row pixel index', fontsize = 14)\n", "ax[0].set_ylabel('Counts', fontsize = 14)\n", "ax[0].set_title('Comparison before/after RefPix step', fontsize = 14)\n", "ax[0].legend()\n", "\n", "ax[1].set_title('Same, but median-substracted counts', fontsize = 14)\n", "ax[1].plot(superbias_results.data[0,-1,:,1500] - np.nanmedian(superbias_results.data[0,-1,:,1500]))\n", "ax[1].plot(refpix_results.data[0,-1,:,1500] - np.nanmedian(refpix_results.data[0,-1,:,1500]))\n", "ax[1].set_xlabel('Row pixel index', fontsize = 14)\n", "ax[1].set_ylabel('Counts - Median Counts', fontsize = 14)" ] }, { "cell_type": "markdown", "id": "b66b2f7d", "metadata": {}, "source": [ "Wow! That looks much better. Several things have been removed, including some interesting high-frequency noise happening in the rows. That's the so-called [odd-even effect](https://github.com/ers-transit/hackaton-2021-jwst-pipeline-interaction/blob/master/Interacting%20with%20the%20JWST%20pipeline%20---%20the%20case%20of%20NIRISS.ipynb), that the `refpix` step takes care of efficiently thanks to reference pixels (pixels insenstive to light) in the detector.\n", "\n", "As can be seen, most of the detector structure is taken care of up to this point. The backgrounds are now nicely suited slightly _above_ zero, as they should; most detector effects are gone and the group looks much cleaner. It is very instructive to do this kind of visual checks on real data, as they can significantly impact the final achieved S/N if not properly accounted for.\n", "\n", "## 3.3 Linearity corrections\n", "\n", "Perhaps one of the most important detector effects to be sure is done well is the linearity correction. This can have dramatic effects on pixels attaining high fluxes; on the order of a couple percent on the precision one measures overall flux changes, which is in the order of the variation of e.g., a transit lightcurve. \n", "\n", "Right away from the uncalibrated data, the up-the-ramp samples in JWST detectors are non-linear, with the pixels at the lower fluences being almost linear and pixels near the saturation ranges deviating significantly from this behavior. For our current dataset, however, this is hard to see directly on a given up-the-ramp sample, though. Take, for instance, the pixel located in row index 45 and column index 1500 (which in this example is one of the pixels that attains the larger signals):" ] }, { "cell_type": "code", "execution_count": null, "id": "4b830de4", "metadata": {}, "outputs": [], "source": [ "# Fit a line to the data:\n", "coeff = np.polyfit(np.arange(6) + 1, refpix_results.data[0, :, 45, 1500], 1)\n", "\n", "plt.figure(figsize = (6,4))\n", "plt.plot(np.arange(6) + 1, refpix_results.data[0, :, 45, 1500], 'o-', label = 'Samples up-the-ramp')\n", "plt.plot(np.arange(6) + 1, np.polyval(coeff, np.arange(6) + 1), 'r--', label = 'Best-fit line')\n", "plt.legend()\n", "plt.xlabel('Group number', fontsize = 14)\n", "plt.ylabel('Counts', fontsize = 14)\n", "plt.title('First integration of pixel (40, 1500)', fontsize = 14)\n", "\n", "plt.figure(figsize = (6,2))\n", "plt.plot(np.arange(6) + 1, refpix_results.data[0, :, 45, 1500] - np.polyval(coeff, np.arange(6) + 1), 'o-')\n", "plt.xlabel('Group number', fontsize = 14)\n", "plt.ylabel('Counts', fontsize = 14)\n", "plt.title('Counts - Best-fit line', fontsize = 14)" ] }, { "cell_type": "markdown", "id": "9670f660", "metadata": {}, "source": [ "For these count levels, and for this particular pixel, the detector seems to behave mostly okay, but this is somewhat of a tricky plot. How much of the non-linearity is hidden in the noise? Could this pixel have just be a special one that is behaving OK? We need a population analysis of all the pixels to figure out if everything looks OK. \n", "\n", "One trick to glance at how the linearity of the up-the-ramp samples evolves as one goes up-the-ramp is to note that if the detector is linear, it doesn't matter at which up-the-ramp sample one looks at, the **fluence level should change from group-to-group at _the same rate_ on average**. So one can quickly investigate if linearity is an issue (and if the pipeline is correctly correcting for it) by:\n", "\n", "1. Taking the difference in fluence between two subsequent groups (say, the last two).\n", "2. Taking the difference in fluence between two _other_ subsequent groups (say, the first two).\n", "3. Take the ratio between those differences.\n", "\n", "If the detector is linear, then all the pixels should fall around a ratio of 1. Do they? Let's try this experiment out. Let's first take the difference of the last two and first two groups for all the pixels of all the integrations --- then take the ratio of those:" ] }, { "cell_type": "code", "execution_count": null, "id": "7a175321", "metadata": {}, "outputs": [], "source": [ "last_pair = refpix_results.data[:, 5, :, :] - refpix_results.data[:, 4, :, :]\n", "first_pair = refpix_results.data[:, 1, :, :] - refpix_results.data[:, 0, :, :]" ] }, { "cell_type": "code", "execution_count": null, "id": "5a060399", "metadata": {}, "outputs": [], "source": [ "ratio = last_pair / first_pair" ] }, { "cell_type": "markdown", "id": "596ec70a", "metadata": {}, "source": [ "Let's now flatten those arrays and plot them as a function of total fluence at the very last group. If linearity weren't an issue, all of these should line around 1:" ] }, { "cell_type": "code", "execution_count": null, "id": "0d14c23c", "metadata": {}, "outputs": [], "source": [ "flattened_ratio = ratio.flatten()\n", "flattened_fluences = refpix_results.data[:, 5, :, :].flatten()" ] }, { "cell_type": "code", "execution_count": null, "id": "e53df40c", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize = (6,4))\n", "plt.plot(flattened_fluences, flattened_ratio, '.', alpha = 0.01, color = 'black')\n", "plt.plot([0,20000], [1., 1.], 'r--')\n", "plt.ylim(0.5,1.5)\n", "plt.xlim(0,20000)\n", "plt.xlabel('Fluence at the last group (counts)', fontsize = 14)\n", "plt.ylabel('(6-5) / (2-1) Group differences', fontsize = 14)\n", "plt.title('No linearity correction', fontsize = 14)" ] }, { "cell_type": "markdown", "id": "5bbcd297", "metadata": {}, "source": [ "Indeed, the data does _not_ line up around 1. So linearity _is_ an issue!\n", "\n", "The next step in the pipeline, the `linearity` step, corrects for this effect. We now correct for it and make sure we use the correct reference file that was used to generate our data:" ] }, { "cell_type": "code", "execution_count": null, "id": "f1c22afd", "metadata": {}, "outputs": [], "source": [ "# Download linearity reference file:\n", "ref_filename = 'jwst_niriss_linearity_0013.fits'\n", "download_link = 'https://stsci.box.com/shared/static/vwkayh6no08kjkxeqktimmvrctambdto.fits'\n", "\n", "# Look if the file has been downloaded. If not, do it:\n", "if os.path.exists(data_folder + ref_filename):\n", " \n", " print(ref_filename + ' already exists. Not downloading.')\n", "\n", "else:\n", " \n", " print('Downloading {}...'.format(ref_filename)+' from '+ download_link)\n", " download_filename = download_file(download_link,\n", " cache=False)\n", " # Rename file:\n", " os.rename(download_filename, data_folder + ref_filename)" ] }, { "cell_type": "code", "execution_count": null, "id": "e9ea1c7c", "metadata": {}, "outputs": [], "source": [ "# Run linearity step:\n", "linearity_results = calwebb_detector1.linearity_step.LinearityStep.call(refpix_results, \n", " output_dir=pipeline_outputs_folder, \n", " override_linearity = data_folder + ref_filename,\n", " save_results=True)" ] }, { "cell_type": "markdown", "id": "47147cb6", "metadata": {}, "source": [ "All right, let's try the same experiment but on the corrected data:" ] }, { "cell_type": "code", "execution_count": null, "id": "0f8002de", "metadata": {}, "outputs": [], "source": [ "corrected_last_pair = linearity_results.data[:, 5, :, :] - linearity_results.data[:, 4, :, :]\n", "corrected_first_pair = linearity_results.data[:, 1, :, :] - linearity_results.data[:, 0, :, :]\n", "corrected_ratio = corrected_last_pair / corrected_first_pair" ] }, { "cell_type": "markdown", "id": "3041b82e", "metadata": {}, "source": [ "Let's plot:" ] }, { "cell_type": "code", "execution_count": null, "id": "08942e0b", "metadata": {}, "outputs": [], "source": [ "flattened_corrected_ratio = corrected_ratio.flatten()\n", "flattened_corrected_fluences = linearity_results.data[:, 5, :, :].flatten()" ] }, { "cell_type": "code", "execution_count": null, "id": "f0f30a7c", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize = (6,4))\n", "plt.plot(flattened_corrected_fluences, flattened_corrected_ratio, '.', alpha = 0.01, color = 'black')\n", "plt.plot([0,20000], [1., 1.], 'r--')\n", "plt.ylim(0.5,1.5)\n", "plt.xlim(0,20000)\n", "plt.xlabel('Fluence at the last group (counts)', fontsize = 14)\n", "plt.ylabel('(6-5) / (2-1) Group differences', fontsize = 14)\n", "plt.title('After linearity correction', fontsize = 14)" ] }, { "cell_type": "markdown", "id": "91161c62", "metadata": {}, "source": [ "That looks **much** better. Plots like this will be fundamental to make with real data, to make sure linearity was properly corrected and it is not an issue with the current dataset. If it is, say, for the very last groups, one could follow a strategy in which one marks those last groups as saturated pixels --- so future pipeline steps don't consider those in the analysis. Better linearity solutions could also be a possible solution on its own for problems like that although most likely solvable on longer time-scales.\n", "\n", "\n", "## 3.4 1/f noise\n", "\n", "One of the last steps before the most computationally expensive steps in the pipeline is the `dark_current` step. This step grabs a reference file that calculates the dark current at each group, and applies the same correction to every integration in the same way. In our experiment here, the impact of performing this step or not is purely cosmetic --- however, in real data, where the spectrum might move sligthly from integration to integration, it would be interesting to test if this step does more harm than good. We'll apply it here nonetheless:" ] }, { "cell_type": "code", "execution_count": null, "id": "d41eade2", "metadata": {}, "outputs": [], "source": [ "# Run superbias step:\n", "darkcurrent_results = calwebb_detector1.dark_current_step.DarkCurrentStep.call(linearity_results, \n", " output_dir=pipeline_outputs_folder, \n", " save_results=True)" ] }, { "cell_type": "markdown", "id": "755efd9f", "metadata": {}, "source": [ "Let's see its impact on products before the dark current correction:" ] }, { "cell_type": "code", "execution_count": null, "id": "47e46eaf", "metadata": {}, "outputs": [], "source": [ "# Plot them:\n", "plt.figure(figsize=(20,5))\n", "im = plt.imshow(linearity_results.data[0,-1,:,:] / np.nanmedian(linearity_results.data[0,-1,:,:]), \\\n", " interpolation = 'None')\n", "im.set_clim(-3,2)\n", "plt.title('Before the DarkCurrent step:')\n", "\n", "# Plot them:\n", "plt.figure(figsize=(20,5))\n", "im = plt.imshow(darkcurrent_results.data[0,-1,:,:] / np.nanmedian(darkcurrent_results.data[0,-1,:,:]), \\\n", " interpolation = 'None')\n", "im.set_clim(-3,2)\n", "plt.title('After the DarkCurrent step:')" ] }, { "cell_type": "markdown", "id": "342f349d", "metadata": {}, "source": [ "This step is the last in the near-infrared detectors that tries to remove any detector-level effects at the group-level, at least in the current version of the pipeline. However, some detector effects still remain --- in particular, the \"infamous\" 1/f noise. This is a type of noise that appears on every group independently, but it _can_ be buried below some cosmetic effects if we were to plot the groups after each other:" ] }, { "cell_type": "code", "execution_count": null, "id": "9294bff5", "metadata": {}, "outputs": [], "source": [ "for i in range(6):\n", " \n", " plt.figure(figsize=(20,5))\n", " im = plt.imshow(darkcurrent_results.data[0,i,:,:], interpolation = 'None')\n", " im.set_clim(-50,50)\n", " plt.title('First integration, group '+str(i+1))" ] }, { "cell_type": "markdown", "id": "94338c17", "metadata": {}, "source": [ "Instead, and to see it more clearly, let's plot the differences between groups --- this should remove the signals which are common to each group, and only leave any residual noise in the background region. \n", "\n", "Let's plot the difference between groups 2 and 1, 3 and 2, and so on:" ] }, { "cell_type": "code", "execution_count": null, "id": "cf1cfcd0", "metadata": {}, "outputs": [], "source": [ "for i in range(5):\n", " \n", " plt.figure(figsize=(20,5))\n", " im = plt.imshow(darkcurrent_results.data[0,i + 1,:,:] - darkcurrent_results.data[0,i,:,:], interpolation = 'None')\n", " im.set_clim(-50,50)\n", " plt.title('First integration, group '+str(i + 2)+' - group '+str(i+1))" ] }, { "cell_type": "markdown", "id": "24430360", "metadata": {}, "source": [ "See those vertical stripes? **That's the so-called \"1/f noise\"**. \n", "\n", "As you can see, this noise is not constant accross groups but, instead, each group samples a \"realization\" of it. The reason why it's called \"1/f noise\" is because what is seen as vertical stripes is, instead, a time-series that occurs _accross_ the pixels in each group (really, in each frame --- but in the case of the readout mode being used here, a frame _is_ a group). \n", "\n", "How this works is that, in reality, to read these up-the-ramp samples the detector uses some electronics that need to read each pixel's signal sequentially. It starts with the pixel in one of the corners, and then moves to the next pixel in the same column taking, in the case of this detector (and most currently supported modes for JWST detectors), 10 microseconds. Then, it jumps to the next in the same column in 10 microseconds, and so on. When it reaches the end of a column after jumping in this case through 256 pixels, the electronics are ordered to wait 120 microseconds exactly before moving to the next column. And then the process repeats. Here's a schematic showing a cartoon of this process:\n", "\n", "![alt text](one-over-f.png \"1/f noise schematic\")\n", "\n", "1/f noise in this context arises because the readout electronics have noise associated with it, and the power-spectral density (PSD) of this noise, if one \"tags\" each pixel with the clocking process described above, has a $1/f^{\\beta}$ shape --- with the power index, $\\beta$, very close to 1. \n", "\n", "We're currently hard at work at STScI characterizing this effect (and thinking/implementing possible ways of account for it), but one easy way of decreasing its impact, at least on white-light lightcurves and/or binned lightcurves is to perform simple column-to-column substraction of non-iluminated pixels. That is --- simply select all non-illuminated pixels in a column, take their mean/median, and remove it for that column. The simplest version of this idea is to just take the median value of each column --- this assumes the number of non-illuminated pixels is much larger than the illuminated ones:" ] }, { "cell_type": "code", "execution_count": null, "id": "743dedc5", "metadata": {}, "outputs": [], "source": [ "def correct_one_over_f(image):\n", " \n", " corrected_image = np.copy(image)\n", " for i in range(image.shape[1]):\n", " \n", " corrected_image[:, i] -= np.nanmedian(corrected_image[:,i])\n", " \n", " return corrected_image " ] }, { "cell_type": "markdown", "id": "9ef884ce", "metadata": {}, "source": [ "Let's try it out on the group differences above:" ] }, { "cell_type": "code", "execution_count": null, "id": "62db2931", "metadata": {}, "outputs": [], "source": [ "for i in range(5):\n", " \n", " plt.figure(figsize=(20,5))\n", " im = plt.imshow(correct_one_over_f(darkcurrent_results.data[0,i + 1,:,:] - darkcurrent_results.data[0,i,:,:]), \n", " interpolation = 'None')\n", " im.set_clim(-50,50)\n", " plt.title('First integration, group '+str(i + 1)+' - group '+str(i) + ' (1/f corrected)')" ] }, { "cell_type": "markdown", "id": "45385f9f", "metadata": {}, "source": [ "This looks **much better, doesn't it**? It is important, however, to not let your eye trick you: this method only partially corrects 1/f noise. In particular, our experiments show it's useful for lowering the impact of this noise at time/length-scales of order $\\sim 10$ columns long or larger only. This means that 1/f-like residual noise is _still_ present in these groups at time/length-scales smaller than that. \n", "\n", "You can check the level of residual 1/f noise directly on your own TSO products by assigning time-stamps to each pixel, and computing the PSD of the resulting time-series either on a given group or for a set of groups. To illustrate how to do this, let's get the PSD of the non-iluminated pixels on each column. Then, let's take the median of all the PSDs of all columns in order to get a representative PSD of the non-iluminated pixels on each column: " ] }, { "cell_type": "code", "execution_count": null, "id": "ce4cec1f", "metadata": {}, "outputs": [], "source": [ "# Set some definitions for the power-spectral densities:\n", "clock_between_pixels = 10 * 1e-6 # 10 microseconds\n", "min_frequency = 1. / (256 * clock_between_pixels) # inverse clock-time of a column\n", "max_frequency = (1. / clock_between_pixels) * 0.5 # same for two pixels -- Nyquist frequency (hence x 0.5)\n", "\n", "# Define the frequency array at which we'll compute the PSDs; set linspace because it's faster:\n", "frequencies = np.linspace(min_frequency, max_frequency, 10000)\n", "\n", "# Compute PSDs on a single group difference. First, correct it (to get same result as above):\n", "corrected_difference = correct_one_over_f(darkcurrent_results.data[0, 1, :, :] - darkcurrent_results.data[0, 0, :, :])\n", "\n", "# Define time-stamps for each pixel in a column:\n", "times = np.arange(256) * clock_between_pixels\n", "\n", "# Now iterate through columns, generating the PSD at each column:\n", "psds = []\n", "for i in range(2048):\n", " \n", " # First, compute median-absolute deviation of the column flux:\n", " column_flux = corrected_difference[:, i]\n", " sigma_mad = 1.48 * np.nanmedian( np.abs( column_flux - np.nanmedian(column_flux) ) )\n", " \n", " # Now, identify all values within 10-sigma from zero --- which will be the un-illuminated pixels:\n", " idx = np.where(np.abs(column_flux) < sigma_mad)\n", " \n", " # Compute the PSD of those pixels:\n", " psds.append(LombScargle(times[idx], column_flux[idx], normalization = 'psd').power(frequencies))\n", " \n", "# Compute median PSDs of all columns:\n", "psd = np.median( np.array(psds), axis = 0)" ] }, { "cell_type": "markdown", "id": "ec769bd9", "metadata": {}, "source": [ "Now let's plot this PSD:" ] }, { "cell_type": "code", "execution_count": null, "id": "ea7bd84e", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize = (10,6))\n", "\n", "# plot PSD in log-space:\n", "ax.plot(frequencies, psd)\n", "ax.set_xscale('log')\n", "ax.set_yscale('log')\n", "ax.set_xlim(400, 40000)\n", "ax.set_ylim(20,4e2)\n", "\n", "# Define properties to define upper axis as well:\n", "f = lambda x: ( 1. / x ) / (10 * 1e-6) \n", "g = lambda x: ( 1. / x ) / (10 * 1e-6)\n", "\n", "ax2 = ax.secondary_xaxis(\"top\", functions=(f, g))\n", "ax2.tick_params(axis='both', which='major', labelsize=13)\n", "ax2.set_xlabel(\"Pixel length-scale\", fontsize = 14)\n", "\n", "# Properties of bottom axis:\n", "ax.set_title('PSD of column background\\n', fontsize = 16)\n", "ax.set_xlabel('Frequency (Hz)', fontsize = 14)\n", "ax.set_ylabel('PSD', fontsize = 14)\n", "ax.tick_params(axis='both', which='major', labelsize=13)" ] }, { "cell_type": "markdown", "id": "576aa487", "metadata": {}, "source": [ "This _definitely_ doesn't look like white-noise (which would be a horizontal line on this plot). Instead, it goes up at lower frequencies with a $1/f$ shape --- hence the name of the effect. \n", "\n", "It is important to realize that 1/f noise is a stochastic process and, thus, it can't really be completely \"corrected\". It can, however, be accounted for in the modelling/spectral extraction in order to account for its known properties. As can be seen from the plot above as well, the smaller the ammount of pixels one considers, the more \"white\" the noise appears to be. For example, it is clear that if one only looks at lengthscales smaller than about 10 pixels, then the underlying noise is mostly white. This means that the smaller the aperture used to extract the spectra, the lower the impact of 1/f noise. A corollary of this is that it also means that the smaller the subarray (and/or the spectral profile), the better for 1/f noise as the noise will look mostly white. \n", "\n", "At the end of the day, it might be that certain apertures will be sweetspots for different datasets. As such, plots like the ones made here are very helfpul to identify how well a given correction for the effect is working and/or to identify limits/properties useful for including in the modelling efforts." ] }, { "cell_type": "markdown", "id": "6c5cafb8", "metadata": {}, "source": [ "## 3.5 The \"jump\" detection\n", "\n", "The next step on the pipeline is the \"jump\" step. What this step tries to detect are \"jumps\" due to cosmic-ray hits in the detector in the up-the-ramp samples. Let's apply it and see its effect on the data directly:" ] }, { "cell_type": "code", "execution_count": null, "id": "c39d9aff", "metadata": {}, "outputs": [], "source": [ "# JWST pipeline outputs are typically of the form initialfile_YourStep.fits. \n", "# Check if jump-step has been already ran:\n", "jump_output_filename = data_folder + filename.split('uncal.fits')[0] + 'jumpstep.fits'\n", "\n", "# If it has, load it. If it hasn't, run it:\n", "if os.path.exists(jump_output_filename):\n", " \n", " print('Jump detection already pre-computed; loading it:')\n", " jump_results = datamodels.RampModel(jump_output_filename)\n", " print('...done!')\n", "\n", "else:\n", "\n", " jump_results = calwebb_detector1.jump_step.JumpStep.call(darkcurrent_results, \n", " output_dir=pipeline_outputs_folder, \n", " save_results=True)" ] }, { "cell_type": "markdown", "id": "206b1c74", "metadata": {}, "source": [ "The position of the detected jumps are identified in the `group_dq` --- any value of 4 means a jump was detected. Let's explore how many of those where detected in the first integration and where:" ] }, { "cell_type": "code", "execution_count": null, "id": "92fa379d", "metadata": {}, "outputs": [], "source": [ "first_integration_jumps = np.where(jump_results.groupdq[0, :, :, :] == 4)" ] }, { "cell_type": "code", "execution_count": null, "id": "fe21647d", "metadata": {}, "outputs": [], "source": [ "print(len(first_integration_jumps[0]), 'pixels with JUMP_DET flag.')" ] }, { "cell_type": "markdown", "id": "5b64bf86", "metadata": {}, "source": [ "Wow, that's a lot of jumps. Let's see _where_ they are located:" ] }, { "cell_type": "code", "execution_count": null, "id": "18ed49d4", "metadata": {}, "outputs": [], "source": [ "# Define location of the jumps:\n", "zero_array = np.zeros( [jump_results.groupdq.shape[2], jump_results.groupdq.shape[3]] )\n", "zero_array[first_integration_jumps[1], first_integration_jumps[2]] = 1.\n", "\n", "# Plot:\n", "plt.figure(figsize=(20,5))\n", "im = plt.imshow(zero_array, interpolation = None)" ] }, { "cell_type": "markdown", "id": "6a487a84", "metadata": {}, "source": [ "Interesting! It is clear the jump detection is, by default, too aggressive: the detected jumps are clearly located in particular regions of the detector. In particular, it is evident the step is detecting outliers following the 1/f banding we observed in the previous subsection. This makes sense with the algorithm the jump detection step uses to detect jumps, which is a two-point difference method, in which [group differences are used to detect jumps](https://jwst-pipeline.readthedocs.io/en/latest/jwst/jump/description.html). Since our group differences are plagged with 1/f noise, then that will impact the jump detection too, which only assumes white read noise and poisson noise.\n", "\n", "There are many ways to deal with this, but the most straightfoward is to simply increase the [threshold for jump detections](https://jwst-pipeline.readthedocs.io/en/latest/jwst/jump/arguments.html). The threshold to define a jump can be easily defined in the pipeline via the `rejection_threshold` parameter --- this is the \"number of sigma\" to call a detection. By default that is 4; let's set it to 10 and try again:" ] }, { "cell_type": "code", "execution_count": null, "id": "8abf04be", "metadata": {}, "outputs": [], "source": [ "jump_output_filename10 = data_folder + filename.split('uncal.fits')[0] + 'threshold_10.fits'\n", "\n", "# If it has, load it. If it hasn't, run it:\n", "if os.path.exists(jump_output_filename10):\n", " \n", " print('Jump detection already pre-computed; loading it:')\n", " jump_results10 = datamodels.RampModel(jump_output_filename10)\n", " print('...done!')\n", "\n", "else:\n", "\n", " # Run jump-detection step:\n", " jump_results10 = calwebb_detector1.jump_step.JumpStep.call(darkcurrent_results, \n", " output_dir=pipeline_outputs_folder, \n", " name = 'threshold_10',\n", " rejection_threshold = 10,\n", " save_results=True) " ] }, { "cell_type": "code", "execution_count": null, "id": "d20e518a", "metadata": {}, "outputs": [], "source": [ "first_integration_jumps = np.where(jump_results10.groupdq[0, :, :, :] == 4)\n", "print(len(first_integration_jumps[0]), 'pixels with JUMP_DET flag.')" ] }, { "cell_type": "code", "execution_count": null, "id": "d823c7b1", "metadata": {}, "outputs": [], "source": [ "# Define location of the jumps:\n", "zero_array = np.zeros( [jump_results10.groupdq.shape[2], jump_results10.groupdq.shape[3]] )\n", "zero_array[first_integration_jumps[1], first_integration_jumps[2]] = 1.\n", "\n", "# Plot:\n", "plt.figure(figsize=(20,5))\n", "im = plt.imshow(zero_array, interpolation = None)" ] }, { "cell_type": "markdown", "id": "e9ea1fad", "metadata": {}, "source": [ "That seems much better! And number of detected jumps has decreased by ~1/3 --- it reaches now about ~3% of pixels. This might not be perfect, and we might want to try a couple of runs of the pipeline end-to-end with different thresholds before deciding on a final one.\n", "\n", "**It is important to note that when actual data comes, this parameter might be one of the most important to play with**. While the algorithm has been tested with simulated data, the final test will come with on-sky data. This will be particularly important to pay attention to during the first few Cycles of JWST science observations: it is possible a single threshold won't be useful for _all_ datasets, and thus checking the results like above will be fundamental. This can have a direct impact on the final signal-to-noise ratio achieved by pipeline data products!\n", "\n", "## 3.6 the `ramp-fitting` step\n", "\n", "The last step of `detector1` is the `ramp_fit` step. This step does something that might _appear_ to be quite simple, but that in reality it's not as trivial as it seems to be: fit a line and get the associated uncertainties to the up-the-ramp samples. Let's go ahead and just do it:" ] }, { "cell_type": "code", "execution_count": null, "id": "86d70ffa", "metadata": {}, "outputs": [], "source": [ "rampfitstep_output0 = data_folder + filename.split('uncal.fits')[0] + 'threshold_10_0_rampfitstep.fits'\n", "rampfitstep_output1 = data_folder + filename.split('uncal.fits')[0] + 'threshold_10_1_rampfitstep.fits'\n", "\n", "if os.path.exists(rampfitstep_output0) and os.path.exists(rampfitstep_output1):\n", " \n", " print('Ramp-fitting already pre-computed; loading it:')\n", " rampfitting_results0 = datamodels.RampModel(rampfitstep_output0)\n", " rampfitting_results1 = datamodels.RampModel(rampfitstep_output1)\n", " rampfitting_results = [rampfitting_results0, rampfitting_results1]\n", " print('Done!')\n", " \n", "else:\n", " rampfitting_results = calwebb_detector1.ramp_fit_step.RampFitStep.call(jump_results10, \n", " output_dir=pipeline_outputs_folder, \n", " save_results=True) " ] }, { "cell_type": "markdown", "id": "96724249", "metadata": {}, "source": [ "All right! Note the products of this step for TSO's is actually a list:" ] }, { "cell_type": "code", "execution_count": null, "id": "b214adf5", "metadata": {}, "outputs": [], "source": [ "print(len(rampfitting_results))" ] }, { "cell_type": "markdown", "id": "159fe870", "metadata": {}, "source": [ "The data associated with the zeroth element of this list (`rampfitting_results[0].data`) has dimensions equal to the size of the frames (rows and columns). The first element (`rampfitting_results[1].data`), has three dimensions, the same as the zeroth but for each integration (we usually refer to this product as the `rateints` product --- i.e., the rates per integration):" ] }, { "cell_type": "code", "execution_count": null, "id": "dc2e9e07", "metadata": {}, "outputs": [], "source": [ "rampfitting_results[0].data.shape" ] }, { "cell_type": "code", "execution_count": null, "id": "990028f3", "metadata": {}, "outputs": [], "source": [ "rampfitting_results[1].data.shape" ] }, { "cell_type": "markdown", "id": "34ae6b9d", "metadata": {}, "source": [ "These results correspond to the \"rates\" (i.e., the slope of the line-fits to the up-the-ramp samples). Let's plot the very first integration to see how the rates for that look like:" ] }, { "cell_type": "code", "execution_count": null, "id": "cf70b1de", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(20,5))\n", "im = plt.imshow(rampfitting_results[1].data[0,:,:])\n", "im.set_clim(-3,3)" ] }, { "cell_type": "markdown", "id": "8998289d", "metadata": {}, "source": [ "That looks OK! However, we still see the strong 1/f banding in these rateint products. What if we perform the simple column-to-column corrections we did at the group-level above, but directly on these `rateint` products?:" ] }, { "cell_type": "code", "execution_count": null, "id": "5d27110f", "metadata": {}, "outputs": [], "source": [ "corrected_rampfit_results = np.copy(rampfitting_results[1].data)\n", "original_rampfit_results = np.copy(rampfitting_results[1].data)\n", "\n", "for integration in range(corrected_rampfit_results.shape[0]):\n", " \n", " corrected_rampfit_results[integration, :, :] = correct_one_over_f( original_rampfit_results[integration, :, :] )\n" ] }, { "cell_type": "code", "execution_count": null, "id": "121ed689", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(20,5))\n", "im = plt.imshow(corrected_rampfit_results[0,:,:])\n", "im.set_clim(-3,3)" ] }, { "cell_type": "markdown", "id": "dc22644d", "metadata": {}, "source": [ "This looks much nicer! However, we do see darker pixels between around columns 1250 and 1800 --- this means there was an overcorrection of the estimated background level. This might very well be the simplistic approach to correct for 1/f: a better approach would be to mask all illuminated pixels, instead of just using the median directly of all pixels in a column. We touch base on this once again on Section 5.2 below. For now, we work with these \"1/f-corrected/reduced\" products." ] }, { "cell_type": "markdown", "id": "65ed9b94", "metadata": {}, "source": [ "4.-A (quick) TSO look at the `Spec2` stage \n", "------------------" ] }, { "cell_type": "markdown", "id": "daf062d4", "metadata": {}, "source": [ "The [second stage of the TSO pipeline](https://jwst-pipeline.readthedocs.io/en/latest/jwst/pipeline/calwebb_spec2.html), following the `detector1`, is the `spec2` stage. In Notebook 1 there was already plenty of discussion about the `spec2` stage and thus, here we will take a different route on the analysis: we will perform the spectral tracing, extraction and lightcurve analysis on our own. This is, by the way, nonetheless needed within the version of the JWST pipeline used in this notebook (`1.3.3`) for NIRISS/SOSS, as there are no good extraction routines implemented for the instrument here.\n", "\n", "Based on the above, we will only use the `assign_wcs` step of the [`spec2`](https://jwst-pipeline.readthedocs.io/en/latest/jwst/pipeline/calwebb_spec2.html) stage. This will allow us to assign various properties to our dataset, including a wavelength map; let's run it and extract this map below:" ] }, { "cell_type": "code", "execution_count": null, "id": "d2f692c6", "metadata": {}, "outputs": [], "source": [ "assign_wcs_results = calwebb_spec2.assign_wcs_step.AssignWcsStep.call(rampfitting_results[1], \n", " output_dir=pipeline_outputs_folder, \n", " save_results=True) " ] }, { "cell_type": "markdown", "id": "791f4632", "metadata": {}, "source": [ "For the sake of time, we will focus here on extracting the wavelength map for order 1. This takes a while so by default we download the resulting map from Box, but if we want to do it directly, the code is also added below. This map will give us the corresponding wavelength to each pixel. One can do this as follows with those `assign_wcs_results`:" ] }, { "cell_type": "code", "execution_count": null, "id": "540e34a8", "metadata": {}, "outputs": [], "source": [ "download_wmap = True" ] }, { "cell_type": "code", "execution_count": null, "id": "bdd68bb7", "metadata": {}, "outputs": [], "source": [ "if download_wmap:\n", " \n", " # Download uncal file:\n", " file_link = 'https://stsci.box.com/shared/static/wijgmg519ksr0nrudgaae7py2hfwd4ud.npy'\n", " wmap_filename = 'wavelength_map.npy'\n", "\n", " # Look if the file has been downloaded. If not, do it:\n", " if os.path.exists(data_folder + wmap_filename):\n", " \n", " print(wmap_filename + ' already exists. Not downloading.')\n", " \n", " wavelength_map = np.load(data_folder + wmap_filename)\n", "\n", " else:\n", " \n", " print('Downloading {}...'.format(wmap_filename))\n", " download_filename = download_file(file_link, cache=False)\n", " # Rename file:\n", " os.rename(download_filename, working_folder + wmap_filename)\n", "\n", " wavelength_map = np.load(working_folder + wmap_filename)\n", " \n", "else: \n", " \n", " order = 1 \n", " rows, columns = assign_wcs_results.data[0,:,:].shape\n", " wavelength_map = np.zeros([rows, columns])\n", "\n", " for row in range(rows):\n", " for column in range(columns):\n", " \n", " wavelength_map[row, column] = assign_wcs_results.meta.wcs(column, row, order)[-1]" ] }, { "cell_type": "markdown", "id": "e0bb723f", "metadata": {}, "source": [ "Let's plot this:" ] }, { "cell_type": "code", "execution_count": null, "id": "b287b73f", "metadata": {}, "outputs": [], "source": [ "fig, ax1 = plt.subplots(figsize=(20, 5))\n", "\n", "im1 = ax1.imshow(wavelength_map, cmap='Reds')\n", "\n", "# Colorbar:\n", "divider = make_axes_locatable(ax1)\n", "cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n", "\n", "cbar = fig.colorbar(im1, shrink = 0.08, cax=cax)\n", "cbar.ax.tick_params(labelsize=13)\n", "cbar.ax.get_yaxis().labelpad = 20\n", "cbar.ax.set_ylabel('Wavelength', rotation=270, fontsize = 19)" ] }, { "cell_type": "markdown", "id": "703c2c9e", "metadata": {}, "source": [ "With this wavelength map and with the data, we can now jump directly to spectral tracing, extraction and lightcurve analyisis:" ] }, { "cell_type": "markdown", "id": "c3ec8cf4", "metadata": {}, "source": [ "5.-General TSO product exploration \n", "------------------\n", "\n", "We now jump to the final part of this analysis, in which we perform some general TSO product exploration. Here, we will first perform some tracing of the spectra; then we will extract the spectra and finally we will do some lightcurve analysis on the segmente of the data under analysis in this notebook. Let's begin!\n", "\n", "\n", "## 5.1 Spectral tracing\n", "\n", "Spectral tracing is the art of figuring out where the spectra _actually_ is. For TSOs, this is particularly important because changes in this position can give rises to flux changes on our time-series due to intra and inter pixel sensitivities.\n", "\n", "We have written a small script to get this done below that works for NIRISS/SOSS data. The premise of the algorithm is that, as we've seen above, NIRISS/SOSS data has an interesting profile which has two peaks; to simplify the tracing algorithm, we transform that two-peaked profile to a single peak through a gaussian filter. Then, we do simple centroiding on each column to figure out where each order actually is. We do this repeatedly and sequentially accross different columns, until we reach the end of the spectrum:" ] }, { "cell_type": "code", "execution_count": null, "id": "8a64d024", "metadata": {}, "outputs": [], "source": [ "def trace_spectrum(image, dqflags, xstart, ystart, profile_radius=20, nsigma = 100, gauss_filter_width=10, xend=None, y_tolerance = 5, verbose = False):\n", " \"\"\"\n", " Function that non-parametrically traces NIRISS/SOSS spectra. First, to get the centroid at xstart and \n", " ystart, it convolves the spatial profile with a gaussian filter, finding its peak through usual flux-weighted \n", " centroiding. Next, this centroid is used as a starting point to find the centroid of the left column through \n", " the same algorithm. \n", " \n", " Parameters\n", " ----------\n", "\n", " image: numpy.array\n", " The image that wants to be traced.\n", " dqflags: ndarray\n", " The data quality flags for each pixel in the image. Only pixels with DQ flags of zero will be used \n", " in the centroiding.\n", " xstart: float\n", " The x-position (column) on which the tracing algorithm will be started\n", " ystart: float\n", " The estimated y-position (row) of the center of the trace. An estimate within 10-20 pixels is enough.\n", " profile_radius: float\n", " Expected radius of the profile measured from its center. Only this region will be used to estimate \n", " the centroids of the spectrum.\n", " nsigma : float\n", " Median filters are applied to each column in search of outliers. This number defines how many n-sigma above the noise level \n", " the residuals of the median filter and the image should be considered outliers.\n", " gauss_filter_width: float\n", " Width of the gaussian filter used to perform the centroiding of the first column\n", " xend: int\n", " x-position at which tracing ends. If none, trace all the columns left to xstart.\n", " y_tolerance: float\n", " When tracing, if the difference between the two difference centroids at two contiguous columns is larger than this, \n", " then assume tracing failed (e.g., cosmic ray).\n", " verbose: boolean\n", " If True, print error messages.\n", "\n", " Returns\n", " -------\n", "\n", " x : numpy.array\n", " Columns at which the centroid is calculated.\n", " y : numpy.array\n", " Calculated centroids.\n", " \"\"\"\n", " \n", " # Define x-axis:\n", " if xend is not None:\n", " x = np.arange(xend, xstart + 1)\n", " else:\n", " x = np.arange(0, xstart + 1)\n", " \n", " # Define y-axis:\n", " y = np.arange(image.shape[0])\n", " \n", " # Define status of good/bad for each centroid:\n", " status = np.full(len(x), True, dtype=bool)\n", " \n", " # Define array that will save centroids at each x:\n", " ycentroids = np.zeros(len(x))\n", " \n", " for i in range(len(x))[::-1]:\n", " xcurrent = x[i]\n", "\n", " # Perform median filter to identify nasty (i.e., cosmic rays) outliers in the column:\n", " mf = median_filter(image[:,xcurrent], size = 5)\n", " residuals = mf - image[:,xcurrent]\n", " mad_sigma = get_mad_sigma(residuals)\n", " column_nsigma = np.abs(residuals) / mad_sigma\n", " \n", " # Extract data-quality flags for current column; index good pixels --- mask nans as well:\n", " idx_good = np.where((dqflags[:, xcurrent] == 0) & (~np.isnan(image[:, xcurrent]) & (column_nsigma < nsigma)))[0] \n", " idx_bad = np.where(~(dqflags[:, xcurrent] == 0) & (~np.isnan(image[:, xcurrent]) & (column_nsigma < nsigma)))[0]\n", " \n", " if len(idx_good) > 0:\n", "\n", " # Replace bad values with the ones in the median filter:\n", " column_data = np.copy(image[:, xcurrent])\n", " column_data[idx_bad] = mf[idx_bad]\n", "\n", " # Convolve column with a gaussian filter; remove median before convolving:\n", " filtered_column = gaussian_filter1d(column_data - \\\n", " np.median(column_data), gauss_filter_width)\n", " \n", " # Find centroid within profile_radius pixels of the initial guess:\n", " idx = np.where(np.abs(y - ystart) < profile_radius)[0]\n", " ycentroids[i] = np.sum(y[idx] * filtered_column[idx]) / np.sum(filtered_column[idx])\n", "\n", " # Get the difference of the current centroid with the previous one (if any):\n", " if xcurrent != x[-1]:\n", "\n", " previous_centroid = ycentroids[i + 1]\n", " difference = np.abs(previous_centroid - ycentroids[i])\n", "\n", " if (difference > y_tolerance):\n", "\n", " if verbose:\n", " print('Tracing failed at column',xcurrent,'; estimated centroid:',ycentroids[i],', previous one:',previous_centroid,'> than tolerance: ',y_tolerance,\\\n", " '. Replacing with closest good trace position.')\n", "\n", " ycentroids[i] = previous_centroid\n", " \n", "\n", " ystart = ycentroids[i]\n", " else:\n", " print(xcurrent,'is a bad column. Setting to previous centroid:')\n", " ycentroids[i] = previous_centroid\n", " status[i] = True\n", " \n", " # Return only good centroids:\n", " idx_output = np.where(status)[0]\n", " return x[idx_output], ycentroids[idx_output]\n", "\n", "def get_mad_sigma(x):\n", "\n", " x_median = np.nanmedian(x)\n", "\n", " return 1.4826 * np.nanmedian( np.abs(x - x_median) )" ] }, { "cell_type": "markdown", "id": "da0e970a", "metadata": {}, "source": [ "As we can see, the code needs an initial guess (`xstart` and `ystart` --- the column and row) on where the spectra is located (and where the tracing starts) and an ending pixel (`xend`) so the tracing stops. The tracing moves \"from right to left\" in our plots above. Let's try tracing Order 1 first.\n", "\n", "To figure those out for this order, let's first try to find `xstart` and `ystart` by plotting the profile of Order 1 just before the reference pixels at the right-end of the spectra on the first integration:" ] }, { "cell_type": "code", "execution_count": null, "id": "a2dddee7", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(10,5))\n", "\n", "plt.title('Profile \"cut\" at column with index 2043')\n", "plt.plot(corrected_rampfit_results[0,:,2043])\n", "plt.xlabel('Row index')\n", "plt.ylabel('Counts')\n", "plt.xlim(50,100)" ] }, { "cell_type": "markdown", "id": "0c9b1e3e", "metadata": {}, "source": [ "All right! So it seems the profile center is around `ystart = 70` if we start at column index `xstart = 2043`. What about `xend`? We want this to be just before the reference pixels at the leftmost part of the profile, so we set `xend = 4`. Let's trace the spectra on all the integrations:" ] }, { "cell_type": "code", "execution_count": null, "id": "586347ac", "metadata": {}, "outputs": [], "source": [ "# This takes tens of seconds, so we load it as well:\n", "\n", "if os.path.exists(data_folder + 'xcentroids1.pkl') and os.path.exists(data_folder + 'ycentroids1.pkl'):\n", " \n", " print('Centroids detected. Extracting...')\n", " x_centroids = pickle.load(open(data_folder + 'xcentroids1.pkl', 'rb'))\n", " y_centroids = pickle.load(open(data_folder + 'ycentroids1.pkl', 'rb'))\n", " print('...done')\n", " \n", "else:\n", "\n", " # We set dictionaries that will save the traced centroids of our spectra:\n", " x_centroids, y_centroids = {}, {}\n", "\n", " # Iterate through all integrations:\n", " for i in range(corrected_rampfit_results.shape[0]):\n", " \n", " x_centroids[i], y_centroids[i] = trace_spectrum(corrected_rampfit_results[i, :, :], \\\n", " np.zeros(corrected_rampfit_results[i, :, :].shape), \\\n", " xstart = 2043, ystart = 70, \\\n", " xend = 4)" ] }, { "cell_type": "markdown", "id": "f6c785e2", "metadata": {}, "source": [ "All right, let's see how these results look like. First, let's plot them directly:" ] }, { "cell_type": "code", "execution_count": null, "id": "51484a39", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(10,5))\n", "for i in range(corrected_rampfit_results.shape[0]):\n", " \n", " plt.plot(x_centroids[i], y_centroids[i], color = 'black', alpha = 0.1)\n", " \n", "plt.xlabel('Column index')\n", "plt.ylabel('Row index')" ] }, { "cell_type": "markdown", "id": "294b13be", "metadata": {}, "source": [ "Those look quite nice! How do they look like on the actual profile?" ] }, { "cell_type": "code", "execution_count": null, "id": "c65886dc", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(20,5))\n", "im = plt.imshow(corrected_rampfit_results[0,:,:])\n", "for i in range(corrected_rampfit_results.shape[0]):\n", " \n", " plt.plot(x_centroids[i], y_centroids[i], color = 'black', alpha = 0.1)\n", " \n", "plt.title('Centroids of Order 1 following spectral trace', fontsize = 14)\n", "im.set_clim(-3,3)" ] }, { "cell_type": "markdown", "id": "e7fab325", "metadata": {}, "source": [ "Perfect! Now, let's pass a smooth function through them, so we can take advantage of the smoothness of the profile to pass through outliers _and_ take out some high-frequency variations in the centroiding. We choose to use a B-spline for this job here and, in particular, select 10 equally spaced knots over the profile. In reality, these knots could be chosen through some cleaver method (e.g., cross-validation), but here we find that this number of knots actually works quite OK:" ] }, { "cell_type": "code", "execution_count": null, "id": "b72e4fee", "metadata": {}, "outputs": [], "source": [ "# First define the number of knots:\n", "nknots = 11 # it's ten in reality because we remove the leftmost edge knot\n", "\n", "# Now we create the array of knots:\n", "xstart, xend = 2043, 4\n", "knots = np.arange(xend, xstart, (xstart - xend) / np.double(nknots))[1:]\n", "\n", "# Add knots at the edges:\n", "#knots = np.append(knots, xstart)\n", "\n", "# And define array over which we'll interpolate with this b-spline:\n", "x = np.arange(xend, xstart + 1)" ] }, { "cell_type": "markdown", "id": "c8ca8cd0", "metadata": {}, "source": [ "Now go through data, interpolate with b-spline, save interpolated values:" ] }, { "cell_type": "code", "execution_count": null, "id": "31abfa86", "metadata": {}, "outputs": [], "source": [ "y_interpolated = {}\n", "\n", "# Iterate through all integrations:\n", "for i in range(corrected_rampfit_results.shape[0]):\n", " \n", " tck = splrep(np.double(x_centroids[i]), y_centroids[i], t = knots)\n", " y_interpolated[i] = splev(x, tck)" ] }, { "cell_type": "markdown", "id": "50afbfce", "metadata": {}, "source": [ "Plot b-splines on top of centroids:" ] }, { "cell_type": "code", "execution_count": null, "id": "dc0132cb", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(15,5))\n", "for i in range(corrected_rampfit_results.shape[0]):\n", " \n", " plt.plot(x_centroids[i], y_centroids[i], color = 'black', alpha = 0.1)\n", " plt.plot(x, y_interpolated[i], '-', color = 'red', alpha = 0.1)\n", "\n", "plt.title('Centroids (black) versus B-splines (red)')\n", "plt.xlabel('Column index')\n", "plt.ylabel('Row index')\n", "plt.xlim(0,2043)" ] }, { "cell_type": "markdown", "id": "c2feb511", "metadata": {}, "source": [ "That looks quite nice! Let's plot them on the profile:" ] }, { "cell_type": "code", "execution_count": null, "id": "222e7da3", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(20,5))\n", "im = plt.imshow(corrected_rampfit_results[0,:,:])\n", "for i in range(corrected_rampfit_results.shape[0]):\n", " \n", " plt.plot(x, y_interpolated[i], color = 'red', alpha = 0.1)\n", "\n", "plt.title(\"B-splines following centroids of Order 1's spectral trace\", fontsize = 14)\n", "im.set_clim(-3,3)" ] }, { "cell_type": "markdown", "id": "7bb9a824", "metadata": {}, "source": [ "**Beautiful!** Let's now repeat this for Order 2. To this end, we actually pass a cutted version of our `rateint` products to the algorithm, in order to not making it get confused with Order 1. We also only trace Order 2 down to pixel 700 --- below this, Order 2's overlap with Order 1 is just too large.\n", "\n", "To figure out where to cut the image, let's plot the profile around column 750, around which it seems one could tell the profiles apart by eye:" ] }, { "cell_type": "code", "execution_count": null, "id": "7d48eac4", "metadata": {}, "outputs": [], "source": [ "plt.plot(corrected_rampfit_results[0,:,750])\n", "plt.xlim(0,120)\n", "plt.title('Profile at column 750')\n", "plt.xlabel('Row index')\n", "plt.ylabel('Counts/sec')" ] }, { "cell_type": "markdown", "id": "c6b2985f", "metadata": {}, "source": [ "So, it seems we could cut the image at row number 70. What abound `xstart` and `ystart` in this case? Let's plot the very last column to the right on which Order 2 signal is still observed:" ] }, { "cell_type": "code", "execution_count": null, "id": "9a3a2d6d", "metadata": {}, "outputs": [], "source": [ "plt.plot(corrected_rampfit_results[0,:,1838])\n", "plt.xlim(0,120)\n", "plt.title('Profile at column 1838')\n", "plt.xlabel('Row index')\n", "plt.ylabel('Counts/sec')\n", "plt.xlim(150,250)\n", "plt.ylim(-10,50)" ] }, { "cell_type": "markdown", "id": "032d6a56", "metadata": {}, "source": [ "It seems if we choose column `xstart = 1838` we should choose `ystart = 220`. If we are passing a cutted version of our products, however, we have to pass `ystart = 220 - cut`, where in our case `cut = 70`. Let's pass our tracing algorithm with these parameters:" ] }, { "cell_type": "code", "execution_count": null, "id": "73a90ddb", "metadata": {}, "outputs": [], "source": [ "# This takes tens of seconds, so we load it as well:\n", "\n", "if os.path.exists(data_folder + 'xcentroids2.pkl') and os.path.exists(data_folder + 'ycentroids2.pkl'):\n", " \n", " print('Centroids detected. Extracting...')\n", " x_centroids2 = pickle.load(open(data_folder + 'xcentroids2.pkl', 'rb'))\n", " y_centroids2 = pickle.load(open(data_folder + 'ycentroids2.pkl', 'rb'))\n", " print('...done')\n", " \n", "else:\n", "\n", " # We set dictionaries that will save the traced centroids of our spectra:\n", " x_centroids2, y_centroids2 = {}, {}\n", "\n", " # Iterate through all integrations:\n", " for i in range(corrected_rampfit_results.shape[0]):\n", " \n", " x_centroids2[i], y_centroids2[i] = trace_spectrum(corrected_rampfit_results[i, 70:, :], \\\n", " np.zeros(corrected_rampfit_results[i, 70:, :].shape), \\\n", " xstart = 1838, ystart = 220 - 70, \\\n", " xend = 700)" ] }, { "cell_type": "markdown", "id": "bb40b2c1", "metadata": {}, "source": [ "All right, let's plot these centroids:" ] }, { "cell_type": "code", "execution_count": null, "id": "187597bb", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(10,5))\n", "for i in range(corrected_rampfit_results.shape[0]):\n", " \n", " plt.plot(x_centroids2[i], y_centroids2[i], color = 'black', alpha = 0.1)\n", " \n", "plt.xlabel('Column index')\n", "plt.ylabel('Row index')" ] }, { "cell_type": "markdown", "id": "70db37d2", "metadata": {}, "source": [ "That looks nice! Let's now plot this on top of the `rateint` products, remembering to add the offset corresponding to the fact that we cutted out a fraction of the products to perform the centroiding:" ] }, { "cell_type": "code", "execution_count": null, "id": "635e8006", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(20,5))\n", "im = plt.imshow(corrected_rampfit_results[0,:,:])\n", "for i in range(corrected_rampfit_results.shape[0]):\n", " \n", " plt.plot(x_centroids2[i], y_centroids2[i] + 70, color = 'black', alpha = 0.1)\n", " \n", "plt.title('Centroids of Order 2 following spectral trace', fontsize = 14)\n", "im.set_clim(-3,3)" ] }, { "cell_type": "markdown", "id": "32120a27", "metadata": {}, "source": [ "That looks good! Let's now ran a B-spline through it:" ] }, { "cell_type": "code", "execution_count": null, "id": "afeca825", "metadata": {}, "outputs": [], "source": [ "# First define the number of knots:\n", "nknots = 11 # it's ten in reality because we remove the leftmost edge knot\n", "\n", "# Now we create the array of knots:\n", "xstart, xend = 1838, 700\n", "knots = np.arange(xend, xstart, (xstart - xend) / np.double(nknots))[1:]\n", "\n", "# Add knots at the edges:\n", "#knots = np.append(knots, xstart)\n", "\n", "# And define array over which we'll interpolate with this b-spline:\n", "x2 = np.arange(xend, xstart + 1)" ] }, { "cell_type": "code", "execution_count": null, "id": "18e6f86e", "metadata": {}, "outputs": [], "source": [ "y_interpolated2 = {}\n", "\n", "# Iterate through all integrations:\n", "for i in range(corrected_rampfit_results.shape[0]):\n", " \n", " tck = splrep(np.double(x_centroids2[i]), y_centroids2[i] + 70, t = knots)\n", " y_interpolated2[i] = splev(x2, tck)" ] }, { "cell_type": "markdown", "id": "df606dcb", "metadata": {}, "source": [ "Let's see how this looks together with the traces for Order 1:" ] }, { "cell_type": "code", "execution_count": null, "id": "3fa06a6e", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(20,5))\n", "im = plt.imshow(corrected_rampfit_results[0,:,:])\n", "for i in range(corrected_rampfit_results.shape[0]):\n", " \n", " plt.plot(x, y_interpolated[i], color = 'red', alpha = 0.1)\n", " plt.plot(x2, y_interpolated2[i], color = 'cornflowerblue', alpha = 0.1)\n", "\n", "plt.title(\"B-splines following centroids of Order 1 and 2 spectral traces\", fontsize = 14)\n", "im.set_clim(-3,3)" ] }, { "cell_type": "markdown", "id": "92a515bc", "metadata": {}, "source": [ "**This looks very nice!**\n", "\n", "Before moving on, one interesting extra detail is that our B-splines actually track _how much the spectral traces shifts as a function of the integration number_ --- for **each** pixel. We can track this through, e.g., the position of the middle pixel as a function of integration number. Let's do this for both orders:" ] }, { "cell_type": "code", "execution_count": null, "id": "8502ad87", "metadata": {}, "outputs": [], "source": [ "# Track middle pixel:\n", "middle_pixel = int(len(x)/2)\n", "middle_pixel2 = int(len(x2)/2)\n", "\n", "# Save position of middle pixel for both orders:\n", "position = np.zeros(corrected_rampfit_results.shape[0])\n", "position2 = np.zeros(corrected_rampfit_results.shape[0])\n", "\n", "for i in range(corrected_rampfit_results.shape[0]):\n", " \n", " position[i] = y_interpolated[i][middle_pixel]\n", " position2[i] = y_interpolated2[i][middle_pixel2]\n", "\n", "plt.figure(figsize=(10,5))\n", "\n", "plt.plot(position - np.median(position), label = 'Order 1', color = 'orangered')\n", "plt.plot(position2 - np.median(position2), label = 'Order 2', color = 'cornflowerblue')\n", "plt.xlabel('Integration index', fontsize = 14)\n", "plt.ylabel('Trace movement (pixel)', fontsize = 14)\n", "plt.title('Trace movement in rows as traced by middle pixel of spectral trace', fontsize = 14)\n", "plt.legend()\n", "plt.xlim(0,len(position))" ] }, { "cell_type": "markdown", "id": "bcde3810", "metadata": {}, "source": [ "In this case, it seems the movement of the trace of both orders is quite small. Indeed, the data was generated without any shifts on the spectral trace --- so this essentially sets the limit to what kind of variations we would be able to see on the trace with JWST data using the methods above. \n", "\n", "Based on our experiment above, it seems we are able to measure pixel position changes at the $\\sim 1/100$th of a pixel for both orders. [Estimates of the impact of sub-pixel drifts show](https://arxiv.org/abs/2010.03576) that drifts of order $\\sim 1/10$th of a pixel could give rise to $\\sim 100$ ppm flux variations that could be corrected through tracking trace movements like this. Scaling these numbers down, tracking position drifts at each pixel, as we've done above, could potentially allow us to track this to correct minute pixel changes down to $\\sim 10$ ppm flux variations. This is a powerful external variable to keep in mind and save for our analyses! \n", "\n", "## 5.2 Coming back to 1/f noise \n", "\n", "Another important applications of tracing the spectra is that now we can mask all the illuminated pixels, and try to perform 1/f corrections like the ones we did in Section 3, but directly masking illuminated pixels thanks to those traces. Let's try to make this correction and see the difference:" ] }, { "cell_type": "code", "execution_count": null, "id": "a91ea241", "metadata": {}, "outputs": [], "source": [ "def correct_one_over_f_with_mask(image, mask):\n", " \"\"\"\n", " Perform 1/f column-to-column substraction but use a mask to define which pixels to use for the median\n", " value of the column. Only use pixels whose value in the mask is 1.\n", " \"\"\"\n", " \n", " corrected_image = np.copy(image)\n", " for i in range(image.shape[1]):\n", " \n", " idx = np.where(mask[:,i] == 1)\n", " corrected_image[:, i] -= np.nanmedian(corrected_image[idx,i])\n", " \n", " return corrected_image " ] }, { "cell_type": "code", "execution_count": null, "id": "2dbdd8bb", "metadata": {}, "outputs": [], "source": [ "corrected_rampfit_results2 = np.copy(rampfitting_results[1].data)\n", "\n", "# All pixels within mask_aperture of the trace get tagged as \"illuminated\":\n", "mask_aperture = 30\n", "\n", "for integration in range(corrected_rampfit_results.shape[0]):\n", " \n", " mask = np.ones( original_rampfit_results[integration, :, :].shape )\n", " \n", " # Mark pixels in Order 1:\n", " for i in range(len(x)):\n", " \n", " column, row = int(x[i]), int(y_interpolated[integration][i])\n", " mask[row - mask_aperture : row + mask_aperture, column] = 0.\n", " \n", " # Mark pixels in Order 2:\n", " for i in range(len(x2)):\n", " \n", " column, row = int(x2[i]), int(y_interpolated2[integration][i])\n", " mask[row - mask_aperture : row + mask_aperture, column] = 0.\n", " \n", " corrected_rampfit_results2[integration, :, :] = correct_one_over_f_with_mask( \n", " original_rampfit_results[integration, :, :],\n", " mask)" ] }, { "cell_type": "code", "execution_count": null, "id": "0ebfbe13", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(20,5))\n", "plt.title('Old method (no trace masking)')\n", "im = plt.imshow(corrected_rampfit_results[0,:,:])\n", "im.set_clim(-3,3)\n", "\n", "plt.figure(figsize=(20,5))\n", "plt.title('New method (with trace masking)')\n", "im = plt.imshow(corrected_rampfit_results2[0,:,:])\n", "im.set_clim(-3,3)" ] }, { "cell_type": "markdown", "id": "20e27d9e", "metadata": {}, "source": [ "Nice! Note how this method in which we expliclty mask illuminated pixels works much better, especially when Order 1 and Order 2 are very separated from each other. In what follows, we continue with those new products; we also leave an **excercise for the reader** to answer with the tools in this notebook:\n", "\n", "**5.2.1. What if we did this 1/f correction at the _group_ level? What would we typically want --- do corrections like these at the rateints level or at the group-level?**\n", "\n", "## 5.3 Spectral extraction \n", "\n", "Next, we perform spectral extraction using our traces. For the sake of time, we will only extract Order 1 --- Order 2 spectral extraction is left as an excercise to the reader. \n", "\n", "To perform this extraction, we use the simple algorithm below, which performs \"simple\" extraction on our traces:" ] }, { "cell_type": "code", "execution_count": null, "id": "e273dbd0", "metadata": {}, "outputs": [], "source": [ "def aperture_extraction(image, x, y, aperture_radius, background_radius=50, error_image=None, correct_bkg=True, method = 'sum'):\n", " \"\"\"\n", " This function takes as inputs two arrays (x,y) that follow the trace,\n", " and returns the added flux over the defined aperture radius (and its error, if an error image\n", " is given as well), substracting in the way any background between the aperture radius and the\n", " background radius. The background is calculated by taking the median of the points between the\n", " aperture_radius and the background_radius.\n", "\n", " Parameters\n", " ----------\n", " image: ndarray\n", " Image from which the spectrum wants to be extracted\n", " x: ndarray\n", " Array with the x-axis of the trace (i.e., the columns, wavelength direction)\n", " y: ndarray\n", " Array with the y-axis of the trace (i.e., rows, spatial direction)\n", " aperture_radius: float\n", " Distance from the center of the trace at which you want to add fluxes.\n", " background_radius: float\n", " Distance from the center of the trace from which you want to calculate the background. The\n", " background region will be between this radius and the aperture_radius.\n", " error_image: ndarray\n", " Image with the errors of each pixel value on the image ndarray above\n", " correct_bkg: boolean\n", " If True, apply background correction. If false, ommit this.\n", " method : string\n", " Method used to perform the extraction. Default is `sum`; `average` takes the average of the non-fractional pixels \n", " used to extract the spectrum. This latter one is useful if the input is a wavelength map.\n", " \"\"\"\n", "\n", " method = method.lower()\n", "\n", " # If average method being used, remove background correction:\n", " if method == 'average':\n", " correct_bkg = False\n", "\n", " # Create array that will save our fluxes:\n", " flux = np.zeros(len(x))\n", " if error_image is not None:\n", " flux_error = np.zeros(len(x))\n", " max_column = image.shape[0] - 1\n", "\n", " for i in range(len(x)):\n", "\n", " # Cut the column with which we'll be working with:\n", " column = image[:,int(x[i])]\n", " if error_image is not None:\n", " variance_column = error_image[:,int(x[i])]**2\n", "\n", " # Define limits given by the aperture_radius and background_radius variables:\n", " if correct_bkg:\n", " left_side_bkg = np.max([y[i] - background_radius, 0])\n", " right_side_bkg = np.min([max_column, y[i] + background_radius])\n", " left_side_ap = np.max([y[i] - aperture_radius, 0])\n", " right_side_ap = np.min([max_column, y[i] + aperture_radius])\n", "\n", " # Extract background, being careful with edges:\n", " if correct_bkg:\n", " bkg_left = column[np.max([0, int(left_side_bkg)]) : np.max([0, int(left_side_ap)])]\n", " bkg_right = column[np.min([int(right_side_ap), max_column]) : np.max([int(right_side_bkg), max_column])]\n", " bkg = np.median(np.append(bkg_left, bkg_right))\n", " else:\n", " bkg = 0.\n", "\n", " # Substract it from the column:\n", " column -= bkg\n", "\n", " # Perform aperture extraction of the background-substracted column, being careful with pixelization\n", " # at the edges. First, deal with left (up) side:\n", " l_decimal, l_integer = np.modf(left_side_ap)\n", " l_integer = int(l_integer)\n", " if l_decimal < 0.5:\n", " l_fraction = (0.5 - l_decimal) * column[np.min([l_integer, max_column])]\n", " l_limit = l_integer + 1\n", " if error_image is not None:\n", " l_fraction_variance = ((0.5 - l_decimal)**2) * variance_column[np.min([l_integer, max_column])]\n", " else:\n", " l_fraction = (1. - (l_decimal - 0.5)) * column[np.min([l_integer + 1, max_column])]\n", " l_limit = l_integer + 2\n", " if error_image is not None:\n", " l_fraction_variance = ((1. - (l_decimal - 0.5))**2) * variance_column[np.min([l_integer + 1, max_column])]\n", "\n", " # Now right (down) side:\n", " r_decimal, r_integer = np.modf(right_side_ap)\n", " r_integer = int(r_integer)\n", " if r_decimal < 0.5:\n", " r_fraction = (1. - (0.5 - r_decimal)) * column[np.min([max_column, r_integer])]\n", " r_limit = r_integer\n", " if error_image is not None:\n", " r_fraction_variance = ((1. - (0.5 - r_decimal))**2) * variance_column[np.min([max_column, r_integer])]\n", " else:\n", " r_fraction = (r_decimal - 0.5) * column[np.min([max_column, r_integer + 1])]\n", " r_limit = r_integer + 1\n", " if error_image is not None:\n", " r_fraction_variance = ((r_decimal - 0.5)**2) * variance_column[np.min([max_column, r_integer + 1])]\n", "\n", " # Save total flux in current column:\n", " if method == 'sum':\n", " flux[i] = l_fraction + r_fraction + np.sum(column[l_limit:r_limit])\n", "\n", " elif method == 'average':\n", " flux[i] = np.mean(column[l_limit:r_limit])\n", "\n", " else:\n", " raise Exception('Method \"'+method+'\" currently not supported for aperture extraction. Select either \"sum\" or \"average\".')\n", "\n", " if error_image is not None:\n", " # Total error is the sum of the variances:\n", " flux_error[i] = np.sqrt(np.sum(variance_column[l_limit:r_limit]) + l_fraction_variance + \\\n", " r_fraction_variance)\n", " if error_image is not None:\n", " return flux, flux_error\n", " else:\n", " return flux" ] }, { "cell_type": "markdown", "id": "62d12186", "metadata": {}, "source": [ "Note the above defined code receives an _error image_ --- i.e., the uncertainties on each of the rates. This is actually given by the pipeline in the same `rateint` product object:" ] }, { "cell_type": "code", "execution_count": null, "id": "8725557a", "metadata": {}, "outputs": [], "source": [ "errors = rampfitting_results[1].err" ] }, { "cell_type": "markdown", "id": "69755ca8", "metadata": {}, "source": [ "Let's use this error \"image\" combined with our 1/f-\"reduced\" products and our traces to get the spectrum for our target. Let's first choose an extraction aperture and \"see\" if we are satisfied with it using the first integration:" ] }, { "cell_type": "code", "execution_count": null, "id": "98172c42", "metadata": {}, "outputs": [], "source": [ "# Define distance from the trace which we'll add to get the spectra:\n", "aperture_radius = 15\n", "\n", "# Plot this aperture:\n", "\n", "plt.figure(figsize=(20,5))\n", "plt.title('Visualization of a '+str(aperture_radius)+'-pixel aperture for extraction')\n", "im = plt.imshow(corrected_rampfit_results2[0,:,:])\n", "\n", "# Order 1:\n", "plt.text(1750, 100, 'Order 1', color = 'white', fontsize = 16)\n", "plt.plot(x, y_interpolated[0], color = 'red')\n", "plt.fill_between(x, y_interpolated[0] - aperture_radius, \\\n", " y_interpolated[0] + aperture_radius, \\\n", " color = 'red', alpha = 0.4)\n", "\n", "# Order 2:\n", "plt.text(1500, 215, 'Order 2', color = 'white', fontsize = 16)\n", "plt.plot(x2, y_interpolated2[0], color = 'cornflowerblue')\n", "plt.fill_between(x2, y_interpolated2[0] - aperture_radius, \\\n", " y_interpolated2[0] + aperture_radius, \\\n", " color = 'cornflowerblue', alpha = 0.4)\n", "\n", "im.set_clim(-3,3)" ] }, { "cell_type": "markdown", "id": "c94c3b36", "metadata": {}, "source": [ "This aperture size seems like a good compromise between little overlap between Order 1 and 2 _and_ containing most of the flux accross the trace. Let's go with it! Let's extract the spectra:" ] }, { "cell_type": "code", "execution_count": null, "id": "f260aa87", "metadata": {}, "outputs": [], "source": [ "spectra = {}\n", "\n", "for integration in range(corrected_rampfit_results2.shape[0]):\n", " \n", " spectra[integration] = {}\n", " \n", " # Extract spectra:\n", " spectra[integration]['flux'], spectra[integration]['flux_error'] = aperture_extraction(\n", " corrected_rampfit_results2[integration, :, :], \\\n", " x, \\\n", " y_interpolated[integration], \\\n", " aperture_radius = aperture_radius, \\\n", " correct_bkg = False, \\\n", " error_image = errors[integration, :, :] )\n", " \n", " # Use wavelength map to extract wavelengths:\n", " spectra[integration]['wavelength'] = aperture_extraction(wavelength_map, \\\n", " x, \\\n", " y_interpolated[integration], \\\n", " aperture_radius = aperture_radius, \\\n", " method = 'average') " ] }, { "cell_type": "markdown", "id": "4eb51eef", "metadata": {}, "source": [ "Let's plot the spectra on top of each other:" ] }, { "cell_type": "code", "execution_count": null, "id": "2539079d", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(10,5))\n", "\n", "for integration in range(corrected_rampfit_results.shape[0]):\n", "\n", " plt.plot(spectra[integration]['wavelength'], spectra[integration]['flux'], color = 'black', alpha = 0.1)\n", "\n", "plt.xlim(np.min(spectra[integration]['wavelength']), np.max(spectra[integration]['wavelength']))\n", "plt.xlabel('Wavelength ($\\mu$m)')\n", "plt.ylabel('Counts/sec')\n", "plt.ylim(0,10000)" ] }, { "cell_type": "markdown", "id": "1a528445", "metadata": {}, "source": [ "Little close-up gives us a surprise:" ] }, { "cell_type": "code", "execution_count": null, "id": "4cca77b0", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(10,5))\n", "\n", "for integration in range(corrected_rampfit_results.shape[0]):\n", "\n", " plt.plot(spectra[integration]['wavelength'], spectra[integration]['flux'], color = 'black', alpha = 0.1)\n", "\n", "plt.xlim(np.min(spectra[integration]['wavelength']), np.max(spectra[integration]['wavelength']))\n", "plt.title('Close-up --- in and out-of-transit spectra')\n", "plt.xlabel('Wavelength ($\\mu$m)')\n", "plt.ylabel('Counts/sec')\n", "plt.ylim(8500,9500)\n", "plt.xlim(1.1,1.2)" ] }, { "cell_type": "markdown", "id": "ab398f9f", "metadata": {}, "source": [ "This looks pretty nice for such a simple algorithmic extraction!\n", "\n", "We do see some outliers here and there, but this is all very normal --- we've added all the flux in the apertures to extract the spectra, _including bad pixels and/or undetected cosmic-rays_. This is why it's important to look at the spectra in this way, to be sure that all is OK. In theory we could go back and try to perform some kind of fancy interpolation to \"fix\" bad pixels --- but we won't do this here, and will move ahead with this.\n", "\n", "## 5.4 Quick lightcurve exploration \n", "\n", "Above, we plotted each of the spectra on top of each other. However, most TSO users will want to take a look at the _flux as a function of time_ at each wavelength, or at a set of wavelenghts (e.g., by adding the flux at several wavelengths forming a \"bin\"). This is typically cumbersome to do; however, a good trick is to create a \"lightcurve image\": a plot in which colors indicate flux variations, and on which the x-axis denotes wavelength, and the y-axis denotes time and/or integration. Let's do a plot like this here.\n", "\n", "First, let's put all the lightcurves at each wavelength in a matrix. We'll normalize the fluxes here by the median fluxes at each wavelength bin --- we also assume the wavelengths don't change accross integrations:" ] }, { "cell_type": "code", "execution_count": null, "id": "f0a5ae1c", "metadata": {}, "outputs": [], "source": [ "matrices = {}\n", "\n", "\n", "matrices['order1 wavelengths'] = spectra[integration]['wavelength']\n", "matrices['order1 lightcurves'] = np.zeros([len(spectra[integration]['wavelength']), \\\n", " corrected_rampfit_results.shape[0]])\n", "matrices['order1 errors'] = np.zeros([len(spectra[integration]['wavelength']), \\\n", " corrected_rampfit_results.shape[0]])\n", "\n", "for wavelength_index in range(len(spectra[integration]['wavelength'])):\n", " \n", " fluxes = np.zeros(corrected_rampfit_results.shape[0])\n", " errors = np.zeros(corrected_rampfit_results.shape[0])\n", " \n", " for integration in range(corrected_rampfit_results.shape[0]):\n", " \n", " fluxes[integration] = spectra[integration]['flux'][wavelength_index]\n", " errors[integration] = spectra[integration]['flux_error'][wavelength_index]\n", "\n", " median_fluxes = np.nanmedian(fluxes)\n", " norm_lc = fluxes / median_fluxes\n", " norm_err = errors / median_fluxes\n", " matrices['order1 lightcurves'][wavelength_index, :] = (norm_lc - 1.)*1e6\n", " matrices['order1 errors'][wavelength_index, :] = norm_err * 1e6" ] }, { "cell_type": "markdown", "id": "1256acbd", "metadata": {}, "source": [ "Let's check the lightcurves all together:" ] }, { "cell_type": "code", "execution_count": null, "id": "05ec7521", "metadata": {}, "outputs": [], "source": [ "for i in range(2040):\n", " plt.plot(matrices['order1 lightcurves'][i, :], color = 'black', alpha = 0.01)\n", "plt.ylim(-30000, 25000)\n", "plt.xlabel('Integration number')\n", "plt.ylabel('Relative flux - 1 (ppm)')" ] }, { "cell_type": "markdown", "id": "4f5f4030", "metadata": {}, "source": [ "That's an egress! Makes sense, as we were taking the second segment of data. The entire event is centered around mid-transit. Let's now see it in full colors:" ] }, { "cell_type": "code", "execution_count": null, "id": "4345c000", "metadata": {}, "outputs": [], "source": [ "fig, ax1 = plt.subplots(figsize=(20, 5))\n", "\n", "# First, order 1:\n", "plt.title('Order 1 lightcurves', fontsize = 20)\n", "plt.xlabel('Wavelength (microns)', fontsize = 19)\n", "plt.ylabel('Integration', fontsize = 19)\n", "im1 = ax1.imshow(matrices['order1 lightcurves'].T, interpolation='none', cmap='Reds')\n", "im1.set_clim(-30000, 20000)\n", "\n", "# Y axis:\n", "ticks = np.arange(0, len(matrices['order1 wavelengths']), 250) \n", "ticklabels = [\"{:0.2f}\".format(matrices['order1 wavelengths'][i]) for i in ticks]\n", "ax1.set_xticks(ticks)\n", "ax1.set_xticklabels(ticklabels, fontsize=13)\n", "\n", "# X axis:\n", "ticks = np.arange(0, corrected_rampfit_results.shape[0], 100) \n", "ticklabels = [\"{:0.2f}\".format(i) for i in ticks]\n", "ax1.set_yticks(ticks)\n", "ax1.set_yticklabels(ticklabels, fontsize=13)\n", "\n", "# Colorbar:\n", "divider = make_axes_locatable(ax1)\n", "cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n", "\n", "cbar = fig.colorbar(im1, shrink = 0.08, cax=cax)\n", "cbar.ax.tick_params(labelsize=13)\n", "cbar.ax.get_yaxis().labelpad = 20\n", "cbar.ax.set_ylabel('ppm', rotation=270, fontsize = 19)" ] }, { "cell_type": "markdown", "id": "95573c38", "metadata": {}, "source": [ "That looks nice! There is clearly an egress event happening around integration $\\sim50$; and it's not obvious there are any bad integrations. Bad pixels impacting the time-series are, however obvious and they occur at particular wavelengths. Note that the wavelengths in this plot go from red to blue --- because they match the wavelengths in the profile (leftmost part of the profile is the red end, blue-end is the rightmost side).\n", "\n", "To see the entire transit event (which was injected to the data), one would need to repeat the analysis in this notebook for the first segment of data. **This is left as an excercise to the reader**:\n", "\n", "**5.4.1. Create the same plot as above but _for the entire dataet_, i.e., reduce the first segment of data and \"paste\" it with this second segment. Are you able to see the entire transit?**\n", "\n", "**5.4.2. Perform the same plot/analysis for Order 2. Are you able to get the transit as well?**\n", "\n", "## 5.5 Quick lightcurve sensitivity exploration \n", "\n", "Another very quick analysis one can make with this TSO data are so-called \"sensitivity\" plots: how sensitive is TSO data in comparison with the calculated errorbars at each wavelength. Experiments like this are particularly useful for datasets on which the duration of the astrophysical event is known (like in the case of this simulated dataset which has an evident transit event), or for which residuals can be obtained after modelling the astrophysical event plus systematics. Then, one can perform analyses like this on the residuals directly.\n", "\n", "We simplify the analysis here and instead use the out-of-event data (out-of-transit data, in our case) and compute the typical scatter of the time-series. Then, we compare that against the average expected errorbars given by the pipeline; in other words:\n", "\n", "$\\textrm{Sensitivity}=\\frac{\\textrm{Scatter of the data}}{\\textrm{Estimated errorbars}}$\n", "\n", "If they match, then the sensitivity is best; we understand the errorbars very well. If they are larger than 1, then we are somehow underestimating the errorbars; overestimation happens when this is smaller than 1. \n", "\n", "Let's perform this for our time-series as a function of wavelength. We consider that integrations 60 and above are all out-of-transit:" ] }, { "cell_type": "code", "execution_count": null, "id": "cceaff82", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(15,5))\n", "\n", "ratios = np.sqrt(np.nanvar(matrices['order1 lightcurves'][:, 60:], axis = 1)) / np.nanmedian(matrices['order1 errors'][:, 60:], axis = 1)\n", "plt.plot(matrices['order1 wavelengths'], ratios, 'o', color ='black')\n", "plt.plot(matrices['order1 wavelengths'], median_filter(ratios, size = 41), '-', \\\n", " lw = 4, color ='orangered', label = 'Median filter')\n", "plt.xlim(np.min(matrices['order1 wavelengths']), np.max(matrices['order1 wavelengths']))\n", "print('Min and max, order 1:',np.min(median_filter(ratios, size = 41)), \\\n", " np.max(median_filter(ratios, size = 41)))\n", "plt.xticks(fontsize=14)\n", "plt.yticks(fontsize=14)\n", "plt.title('Order 1 sensitivity analysis', fontsize = 20)\n", "plt.xlabel('Wavelength (microns)', fontsize = 19)\n", "plt.ylabel(r'$\\times$ Expected per-point scatter', fontsize = 19)\n", "plt.ylim(0.5,3.)\n", "plt.legend()\n" ] }, { "cell_type": "markdown", "id": "3bc3d7e1", "metadata": {}, "source": [ "This plot is telling us that the sensitivity for order 1 is about $\\sim \\times1.4$ to $\\sim \\times 2$ times the expected scatter. This is actually OK for such a simple reduction like the one made here (e.g., no bad pixel correction, no double checks for cosmic rays, etc.), but there is a lot of room for improvement. In particular, note that the sensitivity plot above has two terms --- the intrinsic scatter of the data _and_ the errorbars estimated by the pipeline. In this version of the pipeline, for example, the errorbars as estimated by the pipeline in the ramp-fitting step only assume the error on each group is white-noise plus Poisson. The 1/f noise might render this assumption invalid. \n", "\n", "Getting down to a sensitivity of $\\sim 1$ needs we go back and carefully check what steps could be improved. Some of these improvements are left as **excercies left for the reader below**:\n", "\n", "**5.5.1. Repeat the analysis performing the 1/f correction at the group-level. Do the results change? Why?**\n", "\n", "**5.5.2. Repeat the analysis but now correcting cosmic-rays and/or bad pixels by your method of choice. Does this improve the final sensitivities? Is there an absolute increase in precision?**\n", "\n", "**5.5.3. Try different aperture sizes (smaller, larger); do they significantly change the results?**\n", "\n", "**5.5.4. Try performing the same sensitivity plots but skipping one `detector1` step at a time. What is the step that produces the most discrepancy between the sensitivities/lightcurve shapes?**\n", "\n", "## 5.6 White-light lightcurve residual analyses \n", "\n", "It is unlikely actual JWST data will be as perfect as the one we are dealing with in this excercise. Instrumental systematics --- which have been observed in virtually _all_ observatories in the past --- will be one of the main problems we will have to deal and model. To explore those, a high signal-to-noise lightcurve can give us a very good picture to see at which time-scales those systematics show-up, as well as what instrumental variables (e.g., trace position, focus changes, temperature fluctuations) might be the best to model them. \n", "\n", "One way of generating a high signal-to-noise ratio lightcurve is by summing up ('binning') fluxes from several wavelengths. If using all wavelengths, these are typically dubbed \"white-light\" lightcurves; **let's obtain the \"white-light\" lightcurve for Order 1 by summing up the flux at all wavelengths**. We take care of wavelength bins at which there are deviant pixels (due to, e.g., bad pixels, undetected cosmic-rays, etc.) through interpolation here, but we note that the best is always to do this at the group/rateints level, i.e., perform pixel-by-pixel corrections:" ] }, { "cell_type": "code", "execution_count": null, "id": "54f3b6c0", "metadata": {}, "outputs": [], "source": [ "# First, stack all the spectra in a big matrix --- assume no wavelength shifts accross integrations:\n", "wavelengths = spectra[0]['wavelength']\n", "all_spectra = np.zeros([corrected_rampfit_results.shape[0], len(spectra[0]['flux'])])\n", "all_spectra_err = np.zeros([corrected_rampfit_results.shape[0], len(spectra[0]['flux_error'])])\n", "\n", "# Save all spectra to our big array:\n", "for integration in range(corrected_rampfit_results.shape[0]):\n", "\n", " all_spectra[integration, :] = spectra[integration]['flux']\n", " all_spectra_err[integration, :] = spectra[integration]['flux_error']\n", " \n", "# Now obtain the median spectrum using only out-of-transit spectra:\n", "median_spectrum = np.median(all_spectra[60:, :], axis = 0)\n", "median_spectrum_err = np.median(all_spectra_err[60:, :], axis = 0)" ] }, { "cell_type": "markdown", "id": "7e8508c0", "metadata": {}, "source": [ "To figure out bad pixels in this integrated flux, we create a median filtered version of our median spectrum in the wavelength direction. We do the same for the errors. We also calculate the medians of both of those arrays, so we can then re-scale it when we compare those against each spectra:" ] }, { "cell_type": "code", "execution_count": null, "id": "ce086e38", "metadata": {}, "outputs": [], "source": [ "master_mf = median_filter(median_spectrum, size = 11)\n", "master_mf_err = median_filter(median_spectrum_err, size = 11)" ] }, { "cell_type": "markdown", "id": "3963c57e", "metadata": {}, "source": [ "Now, we use this median filtered spectrum --- we substract it to each individual spectrum, and figure out pixels that deviate by more than 10-sigma from it. We identify those wavelength bins as outliers, and replace those by the value of this median-filtered spectrum:" ] }, { "cell_type": "code", "execution_count": null, "id": "135f3484", "metadata": {}, "outputs": [], "source": [ "corrected_all_spectra = np.zeros([corrected_rampfit_results.shape[0], len(spectra[0]['flux'])])\n", "corrected_all_spectra_err = np.zeros([corrected_rampfit_results.shape[0], len(spectra[0]['flux'])])\n", "\n", "# Save all corrected spectra to our new big array:\n", "for integration in range(corrected_rampfit_results.shape[0]):\n", " \n", " # Get MF of current dataset:\n", " cmf = median_filter(spectra[integration]['flux'], size = 11)\n", " cmf_err = median_filter(spectra[integration]['flux_error'], size = 11)\n", " # Ratio with out-of-transit MF:\n", " ratio = np.median(master_mf / cmf)\n", " ratio_err = np.median(master_mf_err / cmf_err)\n", " # Scale master MF to \"local\" MF using the median of these ratios:\n", " mf = master_mf / ratio\n", " mf_err = master_mf_err / ratio_err\n", " \n", " # Do n-sigma rejection:\n", " nsigma_data = np.abs(spectra[integration]['flux'] - mf) / mf_err\n", " idx = np.where(nsigma_data > 10)[0]\n", " corrected_all_spectra[integration, :] = np.copy(spectra[integration]['flux'])\n", " corrected_all_spectra_err[integration, :] = np.copy(spectra[integration]['flux_error'])\n", " corrected_all_spectra[integration, idx] = mf[idx]\n", " corrected_all_spectra_err[integration, idx] = mf_err[idx]" ] }, { "cell_type": "markdown", "id": "2bdf910b", "metadata": {}, "source": [ "Let's plot this corrected spectrum:" ] }, { "cell_type": "code", "execution_count": null, "id": "77259c2e", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(10,5))\n", "\n", "plt.title('Corrected spectrum')\n", "for integration in range(corrected_rampfit_results.shape[0]):\n", "\n", " plt.plot(wavelengths, corrected_all_spectra[integration, :], color = 'black', alpha = 0.1)\n", "\n", "plt.xlim(np.min(spectra[integration]['wavelength']), np.max(spectra[integration]['wavelength']))\n", "plt.xlabel('Wavelength ($\\mu$m)')\n", "plt.ylabel('Counts/sec')\n", "plt.ylim(0,10000)" ] }, { "cell_type": "markdown", "id": "1bd19fdd", "metadata": {}, "source": [ "That's pretty good! Let's now compute the white-light lightcurves by adding up the fluxes:" ] }, { "cell_type": "code", "execution_count": null, "id": "8957f7c6", "metadata": {}, "outputs": [], "source": [ "wl_lightcurve = np.sum(corrected_all_spectra, axis = 1)\n", "wl_lightcurve_error = np.sqrt(np.sum(corrected_all_spectra_err**2, axis = 1))" ] }, { "cell_type": "markdown", "id": "64a2597d", "metadata": {}, "source": [ "Normalize fluxes:" ] }, { "cell_type": "code", "execution_count": null, "id": "a392215f", "metadata": {}, "outputs": [], "source": [ "out_of_transit_flux = np.median(wl_lightcurve[60:])\n", "wl_lightcurve /= out_of_transit_flux\n", "wl_lightcurve_error /= out_of_transit_flux" ] }, { "cell_type": "code", "execution_count": null, "id": "9b31d9db", "metadata": {}, "outputs": [], "source": [ "plt.errorbar(np.arange(len(wl_lightcurve)), wl_lightcurve, wl_lightcurve_error, fmt = '.')" ] }, { "cell_type": "markdown", "id": "22b7a6d5", "metadata": {}, "source": [ "Let's use the out-of-transit integrations to compute a power-spectral density. First, extract time-stamps from `rateints` products:" ] }, { "cell_type": "code", "execution_count": null, "id": "092ed075", "metadata": {}, "outputs": [], "source": [ "# Extract very last times (corresponding to last segment):\n", "times = uncal_data.int_times['int_mid_BJD_TDB'][-len(wl_lightcurve):]\n", "# Convert to seconds since first integration on this last segment:\n", "times_seconds = (times - times[0]) * 24 * 3600." ] }, { "cell_type": "markdown", "id": "658d056c", "metadata": {}, "source": [ "Now compute PSD of the out-of-transit data:" ] }, { "cell_type": "code", "execution_count": null, "id": "c16c0197", "metadata": {}, "outputs": [], "source": [ "frequencies = np.linspace(1. / (times_seconds[-1] - times_seconds[60]), \\\n", " 0.5 * (1. / np.median(np.diff(times_seconds))), 1000)\n", "\n", "psd = LombScargle(times_seconds[60:], wl_lightcurve[60:], normalization = 'psd').power(frequencies)" ] }, { "cell_type": "code", "execution_count": null, "id": "bf9645ee", "metadata": {}, "outputs": [], "source": [ "plt.plot(frequencies, psd)\n", "plt.xscale('log')\n", "plt.yscale('log')\n", "plt.xlabel('Frequency (Hz)')\n", "plt.ylabel('PSD')" ] }, { "cell_type": "markdown", "id": "56a9c264", "metadata": {}, "source": [ "Apparently nothing too special. There is no specific trends in the PSD that tell us there is any residual noise at a particular space at least at the time-lines we computed the PSD at.\n", "\n", "Another interesting plot to make is to **bin the time-series data at different time-scales and compute how the overall scatter evolves with the bin-size**. If the underlying noise is white, then que should see a $\\sim 1/\\sqrt{N_\\textrm{bin}}$ decrease on the scatter, where $N_\\textrm{bin}$ is the number of binned datapoints. Let's see how this plot looks like in our case:" ] }, { "cell_type": "code", "execution_count": null, "id": "27df6dfd", "metadata": {}, "outputs": [], "source": [ "def bin_data(x, y, n_bin):\n", " \n", " x_bins = []\n", " y_bins = []\n", " y_err_bins = []\n", " \n", " for i in range(0,len(x),n_bin):\n", " \n", " x_bins.append(np.median(x[i:i+n_bin-1]))\n", " y_bins.append(np.median(y[i:i+n_bin-1]))\n", " y_err_bins.append(np.sqrt(np.var(y[i:i+n_bin-1]))/np.sqrt(len(y[i:i+n_bin-1])))\n", " \n", " return np.array(x_bins), np.array(y_bins), np.array(y_err_bins)" ] }, { "cell_type": "code", "execution_count": null, "id": "6a1e3a45", "metadata": {}, "outputs": [], "source": [ "nbins = np.arange(2,30,1)\n", "scatter = np.zeros(len(nbins))\n", "\n", "for i in range(len(nbins)):\n", " \n", " xbin, ybin, ybinerr = bin_data(times_seconds[60:], wl_lightcurve[60:], nbins[i])\n", " scatter[i] = np.sqrt(np.var(ybin))" ] }, { "cell_type": "code", "execution_count": null, "id": "9a544e6d", "metadata": {}, "outputs": [], "source": [ "# Plot:\n", "plt.figure(figsize=(7, 5))\n", "plt.plot(nbins, scatter*1e6, 'o-', color = 'black')\n", "plt.plot(nbins, (np.sqrt(nbins[0]) * 250) / np.sqrt(nbins), '--', label='$Expected 1/\\sqrt{N}$', \\\n", " color = 'cornflowerblue')\n", "plt.legend()\n", "plt.xscale('log')\n", "plt.ylabel('RMS (ppm)')\n", "plt.xlabel('Bin size')\n", "plt.xlim(np.min(nbins),np.max(nbins))" ] }, { "cell_type": "markdown", "id": "b70c6b40", "metadata": {}, "source": [ "The plot, overall, seems to indicate (just as the PSD above) that the noise properties are mostly white indeed. Perhaps there is hint for a stronger _decrease_ of the expected behaviour at larger bins, but this is not too obvious. Overall, it seems our white-light data behaves quite nicely as expected, as our simulated data doesn't have any other effect added to it other than the transit event!\n", "\n", "6.-Final words \n", "------------------\n", "\n", "The analyses performed here on JWST products are, by definition, only a first look at exploring TSO products. This is both triggered by the finite time of this JWebbinar, but also due to the fact that we don't have on-sky data yet to test the entire system for TSO observations. Be aware, thus, that the analyses performed here are only a subset of the possible analyses one can perform, and are in no way the final word on _how_ JWST data _should_ be analyzed. This will be solidified more and more as data comes in from this and future cycles.\n", "\n", "I would like to take the opportunity to thank the entire JWST team behind this notebook, that through testing (with actual detector data, through simulations, etc.) and discussions made this product possible. In particular, to the Time-Series Observations Working Group at STScI, especially to Mike Reagan, Stephan Birkmann, Nikolay Nikolov, Arpita Roy, Loïc Albert and Leonardo Ubeda among others. To the NIRCam IDT members, including Everett Schlawin, Jarron Leisenring, Thomas Greene and Thomas Beatty with which a ton of discussion has been ongoing on different topics touched upon here. To the ERS Transiting Exoplanet team who have provided several venues for discussion and community input. To the JWST team behind the pipeline and the mission itself, including and in no particular order Anton Koekemoer, Alicia Canipe, Jeff Valenti, Karl Gordon, Bryan Hilbert and Joseph Filippazzo. " ] }, { "cell_type": "code", "execution_count": null, "id": "24619c03", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.13" } }, "nbformat": 4, "nbformat_minor": 5 }