diff --git a/notebooks/NIRISS/01_niriss_soss_detector1_reduction.ipynb b/notebooks/NIRISS/01_niriss_soss_detector1_reduction.ipynb index 1600b11..7a66650 100644 --- a/notebooks/NIRISS/01_niriss_soss_detector1_reduction.ipynb +++ b/notebooks/NIRISS/01_niriss_soss_detector1_reduction.ipynb @@ -5,7 +5,7 @@ "id": "3e21d550", "metadata": {}, "source": [ - "\"stsci_logo\" \n" + "\"stsci_logo\"\n" ] }, { @@ -14,10 +14,13 @@ "metadata": {}, "source": [ "\n", + "\n", "# NIRISS/SOSS Notebook 1: Downloading and Calibrating 'uncal' TSO Products\n", - "-----\n", + "\n", + "---\n", "\n", "**Authors**:\n", + "\n", "- **Tyler Baines** | Science Support Analyst | NIRISS Branch | tbaines@stsci.edu\n", "- **Néstor Espinoza** | AURA Assistant Astronomer | Mission Scientist for Exoplanet Science | nespinoza@stsci.edu\n", "- **Aarynn Carter** | AURA Assistant Astronomer | NIRISS Branch | aacarter@stsci.edu\n", @@ -29,27 +32,27 @@ "\n", "\n", "## Table of contents\n", + "\n", "1. [Introduction](#introduction)
\n", - " 1.1 [Purpose of this Notebook](#purpose)
\n", - " 1.2 [Data & Context of the Observations](#data)
\n", + " 1.1 [Purpose of this Notebook](#purpose)
\n", + " 1.2 [Data & Context of the Observations](#data)
\n", "2. [Imports](#imports)
\n", "3. [Downloading & Quick Looks at JWST TSO data](#download)
\n", - " 3.1 [Downloading TSO data from MAST](#mast)
\n", - " 3.2 [Quicklook, pt. I: Target Acquisition](#ta)
\n", - " 3.3 [Quicklook, pt. II: `datamodels` & TSO Science Data Products](#science)
\n", + " 3.1 [Downloading TSO data from MAST](#mast)
\n", + " 3.2 [Quicklook, pt. I: Target Acquisition](#ta)
\n", + " 3.3 [Quicklook, pt. II: `datamodels` & TSO Science Data Products](#science)
\n", "4. [A TSO tour through the `Detector1` stage](#detector1)
\n", - " 4.1 [Checking data quality flags](#dqflags)
\n", - " 4.2 [Identifying saturated pixels](#saturation)
\n", - " 4.3 [Removing detector-level effects: the `superbias` and `refpix` steps](#refpix)
\n", - " 4.4 [Linearity corrections](#linearity)
\n", - " 4.5 [Removing the dark current](#dark-current)
\n", - " 4.6 [Correcting 1/f noise](#one_over_f)
\n", - " 4.7 [Detecting \"jumps\" on up-the-ramp sampling](#jump)
\n", - " 4.8 [Fitting ramps with the `ramp_fit` step](#rampfit)
\n", + " 4.1 [Checking data quality flags](#dqflags)
\n", + " 4.2 [Identifying saturated pixels](#saturation)
\n", + " 4.3 [Removing detector-level effects: the `superbias` and `refpix` steps](#refpix)
\n", + " 4.4 [Linearity corrections](#linearity)
\n", + " 4.5 [Removing the dark current](#dark-current)
\n", + " 4.6 [Correcting 1/f noise](#one_over_f)
\n", + " 4.7 [Detecting \"jumps\" on up-the-ramp sampling](#jump)
\n", + " 4.8 [Fitting ramps with the `ramp_fit` step](#rampfit)
\n", "5. [Final words](#final-words)
\n", "\n", - "\n", - "CRDS Context used: jwst_1225.pmap" + "CRDS Context used: jwst_1225.pmap\n" ] }, { @@ -58,17 +61,18 @@ "metadata": {}, "source": [ "## 1. Introduction \n", + "\n", "
\n", "\n", "### 1.1 Purpose of this Notebook\n", "\n", "In this Notebook, we aim to perform an exploration of Time Series Observations (TSO) products, focusing in particular on products obtained by the [Transiting Exoplanet JWST Early Release Science (ERS) team](https://www.stsci.edu/jwst/science-execution/approved-programs/dd-ers/program-1366) --- a real science dataset that we will reduce starting from the most \"raw\" forms of data products that can be downloaded from MAST. We will learn how to download those products, as well as how to load them and make them interact with the JWST Calibration Pipeline to calibrate them. In a companion Notebook, we then perform spectroscopic analyses on this dataset.\n", "\n", - "### 1.2 Data & Context of the Observations \n", + "### 1.2 Data & Context of the Observations\n", "\n", - "The input data for this Notebook are observations from the ERS [Program 1366](https://www.stsci.edu/jwst/science-execution/program-information) where we will explore observations of the exoplanet WASP-39b obtained with the JWST [NIRISS/SOSS](https://jwst-docs.stsci.edu/jwst-near-infrared-imager-and-slitless-spectrograph/niriss-observing-modes/niriss-single-object-slitless-spectroscopy) mode. This mode uses a unique grism ([GR700XD](https://jwst-docs.stsci.edu/jwst-near-infrared-imager-and-slitless-spectrograph/niriss-instrumentation/niriss-gr700xd-grism)) to seperate the light of a source across three diffraction orders covering a broad wavelength range from 0.6-2.8 $\\mu m$ with a moderate spectral resolution (R $\\approx$ 700 at 1.4 $\\mu m$). \n", + "The input data for this Notebook are observations from the ERS [Program 1366](https://www.stsci.edu/jwst/science-execution/program-information) where we will explore observations of the exoplanet WASP-39b obtained with the JWST [NIRISS/SOSS](https://jwst-docs.stsci.edu/jwst-near-infrared-imager-and-slitless-spectrograph/niriss-observing-modes/niriss-single-object-slitless-spectroscopy) mode. This mode uses a unique grism ([GR700XD](https://jwst-docs.stsci.edu/jwst-near-infrared-imager-and-slitless-spectrograph/niriss-instrumentation/niriss-gr700xd-grism)) to seperate the light of a source across three diffraction orders covering a broad wavelength range from 0.6-2.8 $\\mu m$ with a moderate spectral resolution (R $\\approx$ 700 at 1.4 $\\mu m$).\n", "\n", - "
Notes on the validity of this notebook: It is important to realize that this notebook, as it is, is likely to be quickly outdated as new algorithms and fixes are implemented into the JWST Calibration pipeline, as well as new methodologies and studies update our knowledge of optimally calibrating JWST data products. An up-to-date list of known JWST pipeline issues (some of which we touch on this notebook) can be found on the Known Issues with the JWST Data Products JDox page. In doubt, or for any questions, please contact the JWST Helpdesk!" + "
Notes on the validity of this notebook: It is important to realize that this notebook, as it is, is likely to be quickly outdated as new algorithms and fixes are implemented into the JWST Calibration pipeline, as well as new methodologies and studies update our knowledge of optimally calibrating JWST data products. An up-to-date list of known JWST pipeline issues (some of which we touch on this notebook) can be found on the Known Issues with the JWST Data Products JDox page. In doubt, or for any questions, please contact the JWST Helpdesk!\n" ] }, { @@ -77,10 +81,12 @@ "metadata": {}, "source": [ "## 2. Imports \n", + "\n", "
\n", "\n", "For this demonstration we will need the following packages to be instraalled in your python environemnt:\n", - "1. numpy \n", + "\n", + "1. numpy\n", "2. scipy\n", "3. astropy\n", "4. matplotlib\n", @@ -95,7 +101,8 @@ "conda activate jwst-soss-demo-py3.10\n", "pip install -r requirements-soss.txt\n", "```\n", - "With the working python environment setup lets begin importing the packages" + "\n", + "With the working python environment setup lets begin importing the packages\n" ] }, { @@ -134,7 +141,7 @@ "id": "fb2aa5e9", "metadata": {}, "source": [ - "Lastly, lets configure some of the plotting parameters:" + "Lastly, lets configure some of the plotting parameters:\n" ] }, { @@ -143,12 +150,12 @@ "metadata": {}, "outputs": [], "source": [ - "plt.rcParams['figure.figsize'] = (12, 4)\n", - "plt.rcParams['figure.dpi'] = 100\n", - "plt.rcParams['image.origin'] = 'lower'\n", - "plt.rcParams['image.aspect'] = 'auto'\n", - "plt.rcParams['image.interpolation'] = None\n", - "plt.rcParams['image.cmap'] = 'inferno'" + "plt.rcParams[\"figure.figsize\"] = (12, 4)\n", + "plt.rcParams[\"figure.dpi\"] = 100\n", + "plt.rcParams[\"image.origin\"] = \"lower\"\n", + "plt.rcParams[\"image.aspect\"] = \"auto\"\n", + "plt.rcParams[\"image.interpolation\"] = None\n", + "plt.rcParams[\"image.cmap\"] = \"inferno\"" ] }, { @@ -157,13 +164,14 @@ "metadata": {}, "source": [ "## 3. Downloading JWST TSO data \n", + "\n", "
\n", "\n", "The very first step when it comes to analyzing a JWST dataset is to download that data and perform some quick looks so we know the data quality is acceptable to begin with. Here, we will download the `uncal` products, which are one of the \"raw\"-est forms of dataproducts users can download from MAST. We will perform our data download from MAST using `astroquery.mast` and then use the JWST Calibration pipeline to read and have quicklooks at this data. Let's begin!\n", "\n", - "### 3.1 Downloading TSO data from MAST \n", + "### 3.1 Downloading TSO data from MAST\n", "\n", - "To download JWST data from MAST, we will use the `Observations` function from the `astroquery.mast` library. To do this, we need to indicate the properties of the dataset of interest. For this we need to figure out what instrument, filter, program ID _and_ target was obseved. Options for JWST TSO instruments are `NIRISS/SOSS`, `NIRSPEC/SLIT`, `NIRCAM/GRISM`, `MIRI/SLITLESS`, `NIRSPEC/SLIT`, etc. Here, we search for `NIRISS/SOSS`, and the `CLEAR;GR700XD` filter/grating combination, which corresponds to the dataset we want. We define the proposal ID for the ERS program (`1366`) and the name of the target, `WASP-39`:" + "To download JWST data from MAST, we will use the `Observations` function from the `astroquery.mast` library. To do this, we need to indicate the properties of the dataset of interest. For this we need to figure out what instrument, filter, program ID _and_ target was obseved. Options for JWST TSO instruments are `NIRISS/SOSS`, `NIRSPEC/SLIT`, `NIRCAM/GRISM`, `MIRI/SLITLESS`, `NIRSPEC/SLIT`, etc. Here, we search for `NIRISS/SOSS`, and the `CLEAR;GR700XD` filter/grating combination, which corresponds to the dataset we want. We define the proposal ID for the ERS program (`1366`) and the name of the target, `WASP-39`:\n" ] }, { @@ -175,20 +183,27 @@ "source": [ "# Query MAST data\n", "query_results = Observations.query_criteria(\n", - " instrument_name='NIRISS/SOSS',\n", - " filters='CLEAR;GR700XD',\n", - " proposal_id='1366',\n", - " target_name='WASP-39'\n", + " instrument_name=\"NIRISS/SOSS\",\n", + " filters=\"CLEAR;GR700XD\",\n", + " proposal_id=\"1366\",\n", + " target_name=\"WASP-39\",\n", ")\n", "\n", "# Define columns to display\n", "columns_to_display = [\n", - " 'obs_collection', 'instrument_name', 'filters', 'target_name', 'obs_id',\n", - " 's_ra', 's_dec', 't_exptime', 'proposal_id'\n", + " \"obs_collection\",\n", + " \"instrument_name\",\n", + " \"filters\",\n", + " \"target_name\",\n", + " \"obs_id\",\n", + " \"s_ra\",\n", + " \"s_dec\",\n", + " \"t_exptime\",\n", + " \"proposal_id\",\n", "]\n", "\n", "# Display results\n", - "query_results[columns_to_display].show_in_notebook(display_length=3)\n" + "query_results[columns_to_display].show_in_notebook(display_length=3)" ] }, { @@ -196,7 +211,7 @@ "id": "11f10a82-9cef-4f56-b45d-516ee0a37018", "metadata": {}, "source": [ - "This stores _all_ possible observations in the `query_results` variable. Then, we filter all the products to get only the `SCIENCE`, `UNCAL` data products:" + "This stores _all_ possible observations in the `query_results` variable. Then, we filter all the products to get only the `SCIENCE`, `UNCAL` data products:\n" ] }, { @@ -216,7 +231,7 @@ "id": "16e21ae4-11d2-43c3-b40b-c6484420ccbe", "metadata": {}, "source": [ - "Now we'll filter the results to obtain the uncalibrated data products:" + "Now we'll filter the results to obtain the uncalibrated data products:\n" ] }, { @@ -227,9 +242,7 @@ "outputs": [], "source": [ "uncals = Observations.filter_products(\n", - " data_products, \n", - " productType='SCIENCE',\n", - " productSubGroupDescription='UNCAL'\n", + " data_products, productType=\"SCIENCE\", productSubGroupDescription=\"UNCAL\"\n", ")" ] }, @@ -240,7 +253,7 @@ "metadata": {}, "outputs": [], "source": [ - "uncals[['obs_id', 'productSubGroupDescription', 'size']].show_in_notebook()" + "uncals[[\"obs_id\", \"productSubGroupDescription\", \"size\"]].show_in_notebook()" ] }, { @@ -258,9 +271,9 @@ "id": "b6ffcae5-05f1-440f-a43e-965ce0972196", "metadata": {}, "source": [ - "Note how there are 9 data products. The ones with the lowest `size` values are the Target Aquisition exposures, used to align the telescope with the target of interest. These are _very_ useful to check the quality of the observations (and whether or not they were successful!). The _actual_ TSO data are all the products that follow. Note the latter data are segmented --- this is done in the ground to facilitate the processing of the data. \n", + "Note how there are 9 data products. The ones with the lowest `size` values are the Target Aquisition exposures, used to align the telescope with the target of interest. These are _very_ useful to check the quality of the observations (and whether or not they were successful!). The _actual_ TSO data are all the products that follow. Note the latter data are segmented --- this is done in the ground to facilitate the processing of the data.\n", "\n", - "Let's download all the data, including the Target Acquisition frames, which might be useful to diagnose the quality of observations (this might take some time ~5-7 mins):" + "Let's download all the data, including the Target Acquisition frames, which might be useful to diagnose the quality of observations (this might take some time ~5-7 mins):\n" ] }, { @@ -293,7 +306,7 @@ "\n", "The first set of data products we will have a look at are the Target Acquisition (TA) frames. These are frames that are used to precisely center objects in JWST, so as to correct from any JWST blind pointing errors. These frames are taken before the optical element that disperses the light for our TSO observations is put into place, and before doing any small slews to the actual science targets (which should be in the worst case scenario a few tens of arcseconds away from the science target).\n", "\n", - "Note there are 4 TA frames. [This is expected](https://jwst-docs.stsci.edu/jwst-near-infrared-spectrograph/nirspec-operations/nirspec-target-acquisition/nirspec-wide-aperture-target-acquisition); the usual TA WATA procedure (which is used for TSOs) has one exposure that is used to correct for any pointing errors, and a post-correction TA, which is taken as a \"confirmation\" exposure. We can load those frames with the JWST `datamodels`, which as we will see below are extremely useful models to deal with JWST data, as follows: " + "Note there are 4 TA frames. [This is expected](https://jwst-docs.stsci.edu/jwst-near-infrared-spectrograph/nirspec-operations/nirspec-target-acquisition/nirspec-wide-aperture-target-acquisition); the usual TA WATA procedure (which is used for TSOs) has one exposure that is used to correct for any pointing errors, and a post-correction TA, which is taken as a \"confirmation\" exposure. We can load those frames with the JWST `datamodels`, which as we will see below are extremely useful models to deal with JWST data, as follows:\n" ] }, { @@ -303,10 +316,22 @@ "metadata": {}, "outputs": [], "source": [ - "ta1 = datamodels.RampModel(download_dir + '/mastDownload/JWST/jw01366001001_02101_00001-seg001_nis/jw01366001001_02101_00001-seg001_nis_uncal.fits')\n", - "ta2 = datamodels.RampModel(download_dir + '/mastDownload/JWST/jw01366001001_02101_00002-seg001_nis/jw01366001001_02101_00002-seg001_nis_uncal.fits')\n", - "ta3 = datamodels.RampModel(download_dir + '/mastDownload/JWST/jw01366001001_02101_00003-seg001_nis/jw01366001001_02101_00003-seg001_nis_uncal.fits')\n", - "ta4 = datamodels.RampModel(download_dir + '/mastDownload/JWST/jw01366001001_02101_00004-seg001_nis/jw01366001001_02101_00004-seg001_nis_uncal.fits')\n", + "ta1 = datamodels.RampModel(\n", + " download_dir\n", + " + \"/mastDownload/JWST/jw01366001001_02101_00001-seg001_nis/jw01366001001_02101_00001-seg001_nis_uncal.fits\"\n", + ")\n", + "ta2 = datamodels.RampModel(\n", + " download_dir\n", + " + \"/mastDownload/JWST/jw01366001001_02101_00002-seg001_nis/jw01366001001_02101_00002-seg001_nis_uncal.fits\"\n", + ")\n", + "ta3 = datamodels.RampModel(\n", + " download_dir\n", + " + \"/mastDownload/JWST/jw01366001001_02101_00003-seg001_nis/jw01366001001_02101_00003-seg001_nis_uncal.fits\"\n", + ")\n", + "ta4 = datamodels.RampModel(\n", + " download_dir\n", + " + \"/mastDownload/JWST/jw01366001001_02101_00004-seg001_nis/jw01366001001_02101_00004-seg001_nis_uncal.fits\"\n", + ")\n", "\n", "TAs = [ta1, ta2, ta3, ta4]" ] @@ -316,7 +341,7 @@ "id": "67afe9c7-80e1-427b-a821-e5e62e13687d", "metadata": {}, "source": [ - "The `data` attribute of those `datamodels` (e.g., `ta1.data`) stores the actual data from those frames. Let's check the dimensions of those first to familiarize ourselves with those data products:" + "The `data` attribute of those `datamodels` (e.g., `ta1.data`) stores the actual data from those frames. Let's check the dimensions of those first to familiarize ourselves with those data products:\n" ] }, { @@ -327,7 +352,7 @@ "outputs": [], "source": [ "for ta in TAs:\n", - " print(f'Dimensions of the TA frame of {ta.meta.filename}: {ta.shape}')" + " print(f\"Dimensions of the TA frame of {ta.meta.filename}: {ta.shape}\")" ] }, { @@ -337,7 +362,7 @@ "source": [ "The dimensions come in the form `(integrations, groups, pixel, pixel)` --- so both are 1-integration exposures, of 13 groups each, on a 64x64 pixel frame. This actually is exactly what is expected from WATA TA frames!\n", "\n", - "Let's take a look at the TA frames:" + "Let's take a look at the TA frames:\n" ] }, { @@ -358,21 +383,27 @@ "# Plot the TA images\n", "for i, ta in enumerate(TAs):\n", " axes[0, i].imshow(ta.data[int_index, group_index, :, :], vmin=vmin, vmax=vmax)\n", - " axes[0, i].set_title(f'TA_{i+1}')\n", + " axes[0, i].set_title(f\"TA_{i+1}\")\n", "\n", "# Plot the difference images\n", "for i in range(len(TAs)):\n", " if i == 0:\n", - " diff = TAs[i].data[int_index, group_index, :, :] - TAs[-1].data[int_index, group_index, :, :]\n", + " diff = (\n", + " TAs[i].data[int_index, group_index, :, :]\n", + " - TAs[-1].data[int_index, group_index, :, :]\n", + " )\n", " else:\n", - " diff = TAs[i].data[int_index, group_index, :, :] - TAs[i-1].data[int_index, group_index, :, :]\n", + " diff = (\n", + " TAs[i].data[int_index, group_index, :, :]\n", + " - TAs[i - 1].data[int_index, group_index, :, :]\n", + " )\n", " axes[1, i].imshow(diff, interpolation=None)\n", - " axes[1, i].set_title(f'TA_{i+1} - TA_{i if i > 0 else 4}')\n", + " axes[1, i].set_title(f\"TA_{i+1} - TA_{i if i > 0 else 4}\")\n", "\n", - "# add vertical and horizontal crosshair. \n", + "# add vertical and horizontal crosshair.\n", "for ax in axes.ravel():\n", - " ax.axvline(64//2, ls='--', color='white', lw=0.75)\n", - " ax.axhline(64//2, ls='--', color='white', lw=0.75)\n", + " ax.axvline(64 // 2, ls=\"--\", color=\"white\", lw=0.75)\n", + " ax.axhline(64 // 2, ls=\"--\", color=\"white\", lw=0.75)\n", "\n", "plt.tight_layout()\n", "plt.show()" @@ -383,9 +414,9 @@ "id": "dbae5a66-6fb4-4626-8acd-8e4e86aead7b", "metadata": {}, "source": [ - "Nice! There is a source, although it is rather faint, the difference between TA frames supresses the background variation enable the target source to become more apparent given by the sources positive and negative in the differnece frames. The target is reasonably placed around the center of the frame by the 4th exposure. \n", + "Nice! There is a source, although it is rather faint, the difference between TA frames supresses the background variation enable the target source to become more apparent given by the sources positive and negative in the differnece frames. The target is reasonably placed around the center of the frame by the 4th exposure.\n", "\n", - "We won't need the TA loaded in anymore so lets go ahead and release them from memory. " + "We won't need the TA loaded in anymore so lets go ahead and release them from memory.\n" ] }, { @@ -403,9 +434,9 @@ "id": "6ee3cfe5-009e-4088-a00b-ddf7b0cb212e", "metadata": {}, "source": [ - "### 3.3 Quicklook, pt. II: `datamodels` & TSO Science Data Products \n", + "### 3.3 Quicklook, pt. II: `datamodels` & TSO Science Data Products\n", "\n", - "Next up, we will load the **TSO science** data products so we can interact with them. Once again, we open them through the JWST `datamodels` --- and store all segments of data on lists and explored the dimensions of the data products:" + "Next up, we will load the **TSO science** data products so we can interact with them. Once again, we open them through the JWST `datamodels` --- and store all segments of data on lists and explored the dimensions of the data products:\n" ] }, { @@ -416,10 +447,10 @@ "outputs": [], "source": [ "files = [\n", - " '/mastDownload/JWST/jw01366001001_04101_00001-seg001_nis/jw01366001001_04101_00001-seg001_nis_uncal.fits',\n", - " '/mastDownload/JWST/jw01366001001_04101_00001-seg002_nis/jw01366001001_04101_00001-seg002_nis_uncal.fits',\n", - " '/mastDownload/JWST/jw01366001001_04101_00001-seg003_nis/jw01366001001_04101_00001-seg003_nis_uncal.fits',\n", - " '/mastDownload/JWST/jw01366001001_04101_00001-seg004_nis/jw01366001001_04101_00001-seg004_nis_uncal.fits'\n", + " \"/mastDownload/JWST/jw01366001001_04101_00001-seg001_nis/jw01366001001_04101_00001-seg001_nis_uncal.fits\",\n", + " \"/mastDownload/JWST/jw01366001001_04101_00001-seg002_nis/jw01366001001_04101_00001-seg002_nis_uncal.fits\",\n", + " \"/mastDownload/JWST/jw01366001001_04101_00001-seg003_nis/jw01366001001_04101_00001-seg003_nis_uncal.fits\",\n", + " \"/mastDownload/JWST/jw01366001001_04101_00001-seg004_nis/jw01366001001_04101_00001-seg004_nis_uncal.fits\",\n", "]\n", "\n", "uncal_nis = [datamodels.RampModel(download_dir + file) for file in files]" @@ -433,7 +464,7 @@ "outputs": [], "source": [ "for uncal in uncal_nis:\n", - " print(f'Dimensions of the TSO frame: {uncal}')" + " print(f\"Dimensions of the TSO frame: {uncal}\")" ] }, { @@ -445,7 +476,7 @@ "\n", "
Note on memory usage: Loading data products in lists is very useful, but be aware that in particular for TSOs --- which typically involve large data volumes --- Random-Access Memory (RAM) might be severely impacted. The above loaded data products, for instance, take of order ~5 GB --- and this will only be larger as we run pipeline steps below, which convert data products from, e.g., ints to floats, taking even more space. For a typical TSO, when running the pipeline steps we'll run below, consider on the order of ~50 GB will be used. If you don't have 50GB of RAM they should consider alternatives such as a server, or run files individually
\n", "\n", - "JWST `datamodels` simplify accessing vital information from `fits` files, such as instrument/mode, observation dates, and other header values, without needing to know the exact header keywords. For example, to find the observation date in NIS detector data, the datamodels `search` function quickly locates this information, demonstrating its practicality and ease of use." + "JWST `datamodels` simplify accessing vital information from `fits` files, such as instrument/mode, observation dates, and other header values, without needing to know the exact header keywords. For example, to find the observation date in NIS detector data, the datamodels `search` function quickly locates this information, demonstrating its practicality and ease of use.\n" ] }, { @@ -455,7 +486,7 @@ "metadata": {}, "outputs": [], "source": [ - "uncal_nis[0].search('date')" + "uncal_nis[0].search(\"date\")" ] }, { @@ -463,7 +494,7 @@ "id": "71a034dd-e93e-41e4-a707-cc9e14762e0a", "metadata": {}, "source": [ - "Note how we just added a word similar to what we were looking for, and then this function will take a look to find where similar words are located in the `datamodels` attribute tree. From the above, it seems this information is in `meta.date`. Let's check:" + "Note how we just added a word similar to what we were looking for, and then this function will take a look to find where similar words are located in the `datamodels` attribute tree. From the above, it seems this information is in `meta.date`. Let's check:\n" ] }, { @@ -481,7 +512,7 @@ "id": "01e4eb9b-d3b7-4ae2-8757-c26dbc4a8dd2", "metadata": {}, "source": [ - "It works! Note `date_beg` is the one we would be typically interested in checking (which was when the observations happened). Let's try another one. Suppose we wanted to know the name of the PI of this program. Again, the key word here is `pi` so let's insert that one in the `search` function:" + "It works! Note `date_beg` is the one we would be typically interested in checking (which was when the observations happened). Let's try another one. Suppose we wanted to know the name of the PI of this program. Again, the key word here is `pi` so let's insert that one in the `search` function:\n" ] }, { @@ -491,7 +522,7 @@ "metadata": {}, "outputs": [], "source": [ - "uncal_nis[0].search('pi')" + "uncal_nis[0].search(\"pi\")" ] }, { @@ -499,7 +530,7 @@ "id": "72e8c1d5-d797-443c-a126-8b1a8ba24a9d", "metadata": {}, "source": [ - "So this exists under `meta.program.pi_name`:" + "So this exists under `meta.program.pi_name`:\n" ] }, { @@ -517,7 +548,7 @@ "id": "e21230db-23d8-40cc-8cac-9dcd1c8bf041", "metadata": {}, "source": [ - "As for the science frames, dimensions come in the form `(integrations, groups, pixel, pixel)`. So this is a 158-integration segment, with 9 groups per integration each, on a subarray of dimensions 256 x 2048 --- that sounds about right for NIRISS CLEAR/GR700XD exposures. Let's explore integration number 10 of the last group to get and idea of what NIRISS/SOSS data looks like:" + "As for the science frames, dimensions come in the form `(integrations, groups, pixel, pixel)`. So this is a 158-integration segment, with 9 groups per integration each, on a subarray of dimensions 256 x 2048 --- that sounds about right for NIRISS CLEAR/GR700XD exposures. Let's explore integration number 10 of the last group to get and idea of what NIRISS/SOSS data looks like:\n" ] }, { @@ -537,9 +568,9 @@ "vmax = data.mean() * 2.0\n", "\n", "plt.figure(figsize=(12, 3))\n", - "plt.title('Uncal NIS CLEAR/GR700XD data; first segment, integration 10, last group')\n", + "plt.title(\"Uncal NIS CLEAR/GR700XD data; first segment, integration 10, last group\")\n", "plt.imshow(data, vmin=vmin, vmax=vmax)\n", - "plt.colorbar(label='Counts')\n", + "plt.colorbar(label=\"Counts\")\n", "plt.show()" ] }, @@ -548,7 +579,7 @@ "id": "af6118f3-2ad6-4ab0-a45c-382a6de63c41", "metadata": {}, "source": [ - "The data looks great! The sweeping streaks across the dectector columns is the spectrum of WASP-39 dispersed over 3 spectral orders where the brightest spectrum corresponds to order 1 (primary science) which spans a wavelength range from about 0.8 to 2.8 $\\mu m $, the next brightest corresponds to order 2 (secondary science), and followed by order 3 has the lowest throughput of the three and may not be visible raw images. Some structure in the image is mostly dominated by detector-level effects that will be dealt with in the next section of this notebook. " + "The data looks great! The sweeping streaks across the dectector columns is the spectrum of WASP-39 dispersed over 3 spectral orders where the brightest spectrum corresponds to order 1 (primary science) which spans a wavelength range from about 0.8 to 2.8 $\\mu m $, the next brightest corresponds to order 2 (secondary science), and followed by order 3 has the lowest throughput of the three and may not be visible raw images. Some structure in the image is mostly dominated by detector-level effects that will be dealt with in the next section of this notebook.\n" ] }, { @@ -557,9 +588,10 @@ "metadata": {}, "source": [ "## 4. A TSO tour through the `Detector1` stage \n", + "\n", "
\n", "\n", - "The `uncal` data products we loaded above contain a series of detector systematic effects that we need to remove before our data is ready for science. Now, we will move to calibrating those TSO `uncal` data products, which will take care of most of those effects. \n", + "The `uncal` data products we loaded above contain a series of detector systematic effects that we need to remove before our data is ready for science. Now, we will move to calibrating those TSO `uncal` data products, which will take care of most of those effects.\n", "\n", "To perform this calibration, here we will follow most of the steps outlined in the `calwebb_detector1` or \"Stage 1\" processing described in the [JWST Calibration pipeline documentation](https://jwst-pipeline.readthedocs.io/en/latest/jwst/pipeline/calwebb_detector1.html) --- in particular, the one defined for \"Near-IR\" instruments, such as NIRISS. This Stage 1 processing for Near-IR TSOs is defined by a series of steps, which in order are:\n", "\n", @@ -574,7 +606,7 @@ "9. `ramp_fitting`\n", "10. `gain_scale` (not relevant, there is no impact at all for SOSS data)\n", "\n", - "We will slightly modify and/or add some steps to suit our TSO needs below --- let's get started!" + "We will slightly modify and/or add some steps to suit our TSO needs below --- let's get started!\n" ] }, { @@ -584,11 +616,11 @@ "source": [ "### 4.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. Bad pixels are, in fact, curated by the instrument teams in what we colloquially refer to as a \"bad pixel mask\" --- a mask one can \"attach\" to the data products with the JWST Calibration pipeline. This is exactly what the first step in the pipeline, the Data Quality initialized (`dq_init`) step, does. \n", + "An important component of any TSO analysis is to flag bad pixels, pixels identified as cosmic rays and/or identify saturated pixels. Bad pixels are, in fact, curated by the instrument teams in what we colloquially refer to as a \"bad pixel mask\" --- a mask one can \"attach\" to the data products with the JWST Calibration pipeline. This is exactly what the first step in the pipeline, the Data Quality initialized (`dq_init`) step, does.\n", "\n", "#### 4.1.1 Running & understanding the `dq_init` step\n", "\n", - "Let's run the `dq_init` step on the first segment of our NIS data products: " + "Let's run the `dq_init` step on the first segment of our NIS data products:\n" ] }, { @@ -599,8 +631,10 @@ "outputs": [], "source": [ "# Let's run the DQ init step; first for the first segment of the NIS detector:\n", - "print('Running dq_init on NIS:')\n", - "nis_seg1_dqinit = calwebb_detector1.dq_init_step.DQInitStep.call(uncal_nis[0], save_results=False)" + "print(\"Running dq_init on NIS:\")\n", + "nis_seg1_dqinit = calwebb_detector1.dq_init_step.DQInitStep.call(\n", + " uncal_nis[0], save_results=False\n", + ")" ] }, { @@ -608,7 +642,7 @@ "id": "c4a7c514-6606-4c2d-89eb-e5c60ca96186", "metadata": {}, "source": [ - "All right, data-quality flags have been attached to our uncalibrated data products. To figure out why these are so useful, let's take a look at this bad pixel mask that was attached to our data products; in particular, let's peek at the one attached to the NIS detector products. This mask lives in the `pixeldq` attribute of our products (e.g., `nis_seg1_dqinit.pixeldq`). To familiarize ourselves with this, let's print the dimensions of this mask:" + "All right, data-quality flags have been attached to our uncalibrated data products. To figure out why these are so useful, let's take a look at this bad pixel mask that was attached to our data products; in particular, let's peek at the one attached to the NIS detector products. This mask lives in the `pixeldq` attribute of our products (e.g., `nis_seg1_dqinit.pixeldq`). To familiarize ourselves with this, let's print the dimensions of this mask:\n" ] }, { @@ -626,7 +660,7 @@ "id": "7c75f617-64df-4121-a450-fd2ac27a4a69", "metadata": {}, "source": [ - "As expected, it has the same size as our subarray data. [As per the documentation](https://jwst-pipeline.readthedocs.io/en/latest/jwst/references_general/references_general.html?highlight=data%20quality%20flags#data-quality-flags), most data-quality (DQ) flags should be zero in the subarray; let's plot this array to see how many of them get away from this value and where:" + "As expected, it has the same size as our subarray data. [As per the documentation](https://jwst-pipeline.readthedocs.io/en/latest/jwst/references_general/references_general.html?highlight=data%20quality%20flags#data-quality-flags), most data-quality (DQ) flags should be zero in the subarray; let's plot this array to see how many of them get away from this value and where:\n" ] }, { @@ -637,7 +671,7 @@ "outputs": [], "source": [ "plt.figure(figsize=(12, 3))\n", - "plt.title('Non-zero data-quality values across the subarray')\n", + "plt.title(\"Non-zero data-quality values across the subarray\")\n", "plt.imshow(nis_seg1_dqinit.pixeldq, vmin=-0.5, vmax=0.5)\n", "plt.show()" ] @@ -647,7 +681,7 @@ "id": "0d7cbb8f-3281-4bb5-aee2-d1d85bed82e4", "metadata": {}, "source": [ - "All right --- so there are \"special\" pixels all over the place! But, what are the `pixeldq` values telling us? Let's print the pixel in the very corner of the subarray:" + "All right --- so there are \"special\" pixels all over the place! But, what are the `pixeldq` values telling us? Let's print the pixel in the very corner of the subarray:\n" ] }, { @@ -669,11 +703,11 @@ "\n", "#### 4.1.2 Dynamically translating data-quality flags to human-readable form\n", "\n", - "Looking back and forth from the documentation page the data-quality flag values we read from our data-products is a very tedious task. In addition, as we will see below, a pixel can have eventually several flags (e.g., saturated, has a cosmic-ray, etc.) which will, in turn, change some of its data-quality flags to account for this. \n", + "Looking back and forth from the documentation page the data-quality flag values we read from our data-products is a very tedious task. In addition, as we will see below, a pixel can have eventually several flags (e.g., saturated, has a cosmic-ray, etc.) which will, in turn, change some of its data-quality flags to account for this.\n", "\n", "A handy function to convert those data-quality flag numbers to \"human-readable\" form is actually inside the `datamodels` class --- the `datamodels.dqflags`. This simply takes in a data-quality value, and spits out a `set` with strings defining what this is telling us given a so-called \"mnemonic map\" --- one which is actually already loaded in the `datamodels.dqflags.pixel` dictionary.\n", "\n", - "Let's try it out on the data-quality value we observed above:" + "Let's try it out on the data-quality value we observed above:\n" ] }, { @@ -683,7 +717,9 @@ "metadata": {}, "outputs": [], "source": [ - "datamodels.dqflags.dqflags_to_mnemonics(2147483648, mnemonic_map=datamodels.dqflags.pixel)" + "datamodels.dqflags.dqflags_to_mnemonics(\n", + " 2147483648, mnemonic_map=datamodels.dqflags.pixel\n", + ")" ] }, { @@ -691,7 +727,7 @@ "id": "38e09913-917e-4be6-ac14-384a0545696a", "metadata": {}, "source": [ - "Indeed, we get back what we knew --- that is a reference pixel! With this handy-dandy function, we can write a simple snippet to figure out the total tally of all bad pixels as follows:" + "Indeed, we get back what we knew --- that is a reference pixel! With this handy-dandy function, we can write a simple snippet to figure out the total tally of all bad pixels as follows:\n" ] }, { @@ -708,13 +744,12 @@ "\n", "# Iterate through every row and column:\n", "for row in range(rows):\n", - " \n", + "\n", " for column in range(columns):\n", "\n", " # Extract the bad pixel flag(s) for the current pixel at (row, column):\n", " bps = datamodels.dqflags.dqflags_to_mnemonics(\n", - " nis_seg1_dqinit.pixeldq[row, column], \n", - " mnemonic_map=datamodels.dqflags.pixel\n", + " nis_seg1_dqinit.pixeldq[row, column], mnemonic_map=datamodels.dqflags.pixel\n", " )\n", "\n", " # Iterate through the possible flags (it can be more than one!):\n", @@ -752,7 +787,7 @@ "print(separator)\n", "for bp in bad_pixels.keys():\n", " val = bad_pixels[bp]\n", - " frac = 100*(val/float(total_pixels))\n", + " frac = 100 * (val / float(total_pixels))\n", " print(f\"| {bp:<20} | {val:^19} | {frac:^25.2f} |\")\n", "print(separator)" ] @@ -762,9 +797,9 @@ "id": "9e23c07a", "metadata": {}, "source": [ - "Based on our discussion above, we can see some number above make sense. For instance, 10208 `REFERENCE_PIXELS` makes sense as there are a total of 8 columns (4 columns to left and right, and 4 rows to the top of the subarray) --- given the subarray height is 256 pixels and width 2048 pixels (2040 pixels to avoid double counting pixel in the upper corners) given a total of 10208 reference pixels as expected. \n", + "Based on our discussion above, we can see some number above make sense. For instance, 10208 `REFERENCE_PIXELS` makes sense as there are a total of 8 columns (4 columns to left and right, and 4 rows to the top of the subarray) --- given the subarray height is 256 pixels and width 2048 pixels (2040 pixels to avoid double counting pixel in the upper corners) given a total of 10208 reference pixels as expected.\n", "\n", - "Let's go ahead now and attach this bad pixel mask to all the segments of data:" + "Let's go ahead now and attach this bad pixel mask to all the segments of data:\n" ] }, { @@ -784,7 +819,7 @@ "id": "5f5516bd-43b5-4e0a-8223-3a4a31d633ff", "metadata": {}, "source": [ - "
Note on saving data products with the JWST Calibration Pipeline: Sometimes, one might find it useful to save data products after running each step into .fits files, so we can have \"intermediate steps\" stored in our system that we can check at a later time. This can be done when running any of the steps by adding the save_results = True flag to the step calls, e.g., calwebb_detector1.dq_init_step.DQInitStep.call(uncal_nis[i], save_results = True). An output directory can also be defined by using the output_dir parameter." + "
Note on saving data products with the JWST Calibration Pipeline: Sometimes, one might find it useful to save data products after running each step into .fits files, so we can have \"intermediate steps\" stored in our system that we can check at a later time. This can be done when running any of the steps by adding the save_results = True flag to the step calls, e.g., calwebb_detector1.dq_init_step.DQInitStep.call(uncal_nis[i], save_results = True). An output directory can also be defined by using the output_dir parameter.\n" ] }, { @@ -798,7 +833,7 @@ "\n", "#### 4.2.1 Running and understanding the `saturation` step\n", "\n", - "Through the analysis of calibration datasets, the JWST instrument teams have defined signal values for each pixel above which they are considered as \"saturated\". This identification is done through the `saturation` step --- the next step of the JWST pipeline for Detector 1. Let's run it for the very first segment of data for NIS:" + "Through the analysis of calibration datasets, the JWST instrument teams have defined signal values for each pixel above which they are considered as \"saturated\". This identification is done through the `saturation` step --- the next step of the JWST pipeline for Detector 1. Let's run it for the very first segment of data for NIS:\n" ] }, { @@ -817,13 +852,13 @@ "id": "a75bfef9", "metadata": {}, "source": [ - "The saturation step works by primarily comparing the observed count 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_0015.fits` file --- but this might change as new analyses are made and the reference files get updated. \n", + "The saturation step works by primarily comparing the observed count 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_0015.fits` file --- but this might change as new analyses are made and the reference files get updated.\n", "\n", "In addition, at the time of writing, the `saturation` step in the JWST Calibration pipeline [by default flags not only pixels that exceed the signal limit defined by the instrument teams but also all `n_pix_grow_sat` pixels around it](https://jwst-docs.stsci.edu/jwst-science-calibration-pipeline-overview/jwst-operations-pipeline-build-information/jwst-operations-pipeline-build-8-0-release-notes#JWSTOperationsPipelineBuild8.0ReleaseNotes-charge_spilling); which at the time of writing is set to a default of `1`. That means that if a given pixel exceeds the signal limit, all 8 pixels around it will be marked as saturated as well. This is done because it has been observed that \"charge spilling\" can be an issue --- i.e., charge going from one pixel to another. While such migration of charge happens at a wide range of count levels, this is particularly dramatic when a pixel saturates --- reason by which this is set in the pipeline.\n", "\n", "We can check which pixels are saturated in a similar way as to how we checked the data-quality flags in [Section 3.1](#dqflags). The only difference with that analysis is that saturated pixels are integration and group-dependant, i.e., a property of a given pixel _in a given integration and group_. In other words, a pixel that is saturated in one integration and group might have \"recovered\" by the next integration and group.\n", "\n", - "To figure out the data-quality for all integrations and all groups we look at the `groupdq` attribute of our data products instead of the `pixeldq` which we used above. To familiarize ourselves with this, let's print the dimensions of this array first:" + "To figure out the data-quality for all integrations and all groups we look at the `groupdq` attribute of our data products instead of the `pixeldq` which we used above. To familiarize ourselves with this, let's print the dimensions of this array first:\n" ] }, { @@ -841,11 +876,11 @@ "id": "f3cca152-fdd6-4f57-97df-a770da5eb613", "metadata": {}, "source": [ - "As expected, it has dimensions `(integrations, groups, row pixels, column pixels)`, just like the `data` array. The flags in the `groupdq` array follow the same structure as [all the data-quality flags described in the documentation](https://jwst-pipeline.readthedocs.io/en/latest/jwst/references_general/references_general.html?highlight=data%20quality%20flags#data-quality-flags). \n", + "As expected, it has dimensions `(integrations, groups, row pixels, column pixels)`, just like the `data` array. The flags in the `groupdq` array follow the same structure as [all the data-quality flags described in the documentation](https://jwst-pipeline.readthedocs.io/en/latest/jwst/references_general/references_general.html?highlight=data%20quality%20flags#data-quality-flags).\n", "\n", "#### 4.2.2 Exploring saturated pixels via the `groupdq` array\n", "\n", - "To illustrate how to use the `groupdq`, let's pick the last group of integration 10 again and see if any pixels seem to be saturated --- we also count all of the saturated pixels:" + "To illustrate how to use the `groupdq`, let's pick the last group of integration 10 again and see if any pixels seem to be saturated --- we also count all of the saturated pixels:\n" ] }, { @@ -865,21 +900,25 @@ "\n", "verbose = False\n", "for row in range(rows):\n", - " \n", + "\n", " for column in range(columns):\n", "\n", " # Extract the bad pixel flag(s) for the current pixel at (row, column):\n", " bps = datamodels.dqflags.dqflags_to_mnemonics(\n", - " saturation_results.groupdq[integration, group, row, column], \n", - " mnemonic_map=datamodels.dqflags.pixel\n", + " saturation_results.groupdq[integration, group, row, column],\n", + " mnemonic_map=datamodels.dqflags.pixel,\n", " )\n", - " \n", + "\n", " # Check if pixel is saturated; if it is...\n", - " if 'SATURATED' in bps:\n", + " if \"SATURATED\" in bps:\n", "\n", " # ...print which pixel it is, and...\n", " if verbose:\n", - " print('Pixel ({0:},{1:}) is saturated in integration 10, last group'.format(row, column))\n", + " print(\n", + " \"Pixel ({0:},{1:}) is saturated in integration 10, last group\".format(\n", + " row, column\n", + " )\n", + " )\n", "\n", " # ...count it:\n", " nsaturated += 1\n", @@ -887,9 +926,9 @@ " column_idx.append(column)\n", " row_idx.append(row)\n", "\n", - "print(f\"\\nA total of {nsaturated} out of {rows*columns} pixels ({100 * nsaturated / float(rows * columns):.2f}%) are saturated\")\n", - "\n", - "\n" + "print(\n", + " f\"\\nA total of {nsaturated} out of {rows*columns} pixels ({100 * nsaturated / float(rows * columns):.2f}%) are saturated\"\n", + ")" ] }, { @@ -897,7 +936,7 @@ "id": "5b0a4728", "metadata": {}, "source": [ - "As can be seen, not many pixels are saturated on a given group. Let's see how the up-the-ramp samples look like for one of those pixels --- let's say, pixel `(176, 1503)`. Let's show in the same plot the group data-quality flags at each group:" + "As can be seen, not many pixels are saturated on a given group. Let's see how the up-the-ramp samples look like for one of those pixels --- let's say, pixel `(176, 1503)`. Let's show in the same plot the group data-quality flags at each group:\n" ] }, { @@ -911,25 +950,29 @@ "pixel_row, pixel_column = 176, 1503\n", "\n", "plt.figure(figsize=(7, 4))\n", - "plt.title(f'Saturated Pixel: ({pixel_row}, {pixel_column})')\n", - "plt.plot(np.arange(saturation_results.data.shape[1])+1, \n", - " saturation_results.data[integration, :, pixel_row, pixel_column], \n", - " 'o-', color='tomato'\n", + "plt.title(f\"Saturated Pixel: ({pixel_row}, {pixel_column})\")\n", + "plt.plot(\n", + " np.arange(saturation_results.data.shape[1]) + 1,\n", + " saturation_results.data[integration, :, pixel_row, pixel_column],\n", + " \"o-\",\n", + " color=\"tomato\",\n", ")\n", "\n", - "plt.xlim(0.5, saturation_results.data.shape[1]+1.5)\n", - "plt.xlabel('Group number', fontsize=16)\n", - "plt.ylabel('Counts', fontsize=16, color='tomato')\n", + "plt.xlim(0.5, saturation_results.data.shape[1] + 1.5)\n", + "plt.xlabel(\"Group number\", fontsize=16)\n", + "plt.ylabel(\"Counts\", fontsize=16, color=\"tomato\")\n", "\n", "plt.twinx()\n", "\n", - "plt.plot(np.arange(saturation_results.data.shape[1])+1, \n", - " saturation_results.groupdq[integration, :, pixel_row, pixel_column], \n", - " 'o-', color='cornflowerblue'\n", + "plt.plot(\n", + " np.arange(saturation_results.data.shape[1]) + 1,\n", + " saturation_results.groupdq[integration, :, pixel_row, pixel_column],\n", + " \"o-\",\n", + " color=\"cornflowerblue\",\n", ")\n", "\n", - "plt.xlim(0.5, saturation_results.data.shape[1]+1.5)\n", - "plt.ylabel('Group Data-quality', fontsize=16, color='cornflowerblue')\n", + "plt.xlim(0.5, saturation_results.data.shape[1] + 1.5)\n", + "plt.ylabel(\"Group Data-quality\", fontsize=16, color=\"cornflowerblue\")\n", "\n", "plt.show()" ] @@ -939,7 +982,7 @@ "id": "994c506c-e069-4f37-89da-315b8d1a6231", "metadata": {}, "source": [ - "Very interesting plot! Note that all groups appear to be saturated after group ~6 in this example. Likely a cosmic-ray hit happened at this group which left the pixel at a very high count number from group 6 up to the end of the ramp." + "Very interesting plot! Note that all groups appear to be saturated after group ~6 in this example. Likely a cosmic-ray hit happened at this group which left the pixel at a very high count number from group 6 up to the end of the ramp.\n" ] }, { @@ -949,11 +992,11 @@ "source": [ "#### 4.2.3 Setting custom saturation limits with the `saturation` reference file\n", "\n", - "TSOs often obtain data from bright stars that might quickly (i.e., first few groups) give rise to saturated pixels. As described in some early JWST results (see, e.g., [Rustamkulov et al., 2023](https://www.nature.com/articles/s41586-022-05677-y)), in some cases one might even want to be a bit more aggressive on the level of saturation allowed in a given dataset in order to improve on the reliability of the results. As such, understanding how to modify the level of saturation allowed in a given dataset might turn out to be an important skill on real TSO data analysis. \n", + "TSOs often obtain data from bright stars that might quickly (i.e., first few groups) give rise to saturated pixels. As described in some early JWST results (see, e.g., [Rustamkulov et al., 2023](https://www.nature.com/articles/s41586-022-05677-y)), in some cases one might even want to be a bit more aggressive on the level of saturation allowed in a given dataset in order to improve on the reliability of the results. As such, understanding how to modify the level of saturation allowed in a given dataset might turn out to be an important skill on real TSO data analysis.\n", "\n", - "The key file that sets the limits used to call a pixel \"saturated\" is the reference file of the `saturation` step. \n", + "The key file that sets the limits used to call a pixel \"saturated\" is the reference file of the `saturation` step.\n", "\n", - "As discussed above, this can be seen directly on the outputs of the `saturation` step while its running, but it's also saved in our data products:" + "As discussed above, this can be seen directly on the outputs of the `saturation` step while its running, but it's also saved in our data products:\n" ] }, { @@ -971,7 +1014,7 @@ "id": "4e277e99", "metadata": {}, "source": [ - "We can actually load this reference file using the `SaturationModel` as follows:" + "We can actually load this reference file using the `SaturationModel` as follows:\n" ] }, { @@ -982,10 +1025,12 @@ "outputs": [], "source": [ "# Base directory where reference files are stored (this was defined in the Setup section above):\n", - "base_ref_files = os.environ[\"CRDS_PATH\"]+\"/references/jwst/niriss/\"\n", + "base_ref_files = os.environ[\"CRDS_PATH\"] + \"/references/jwst/niriss/\"\n", "\n", "# Read it in:\n", - "saturation_ref_file = datamodels.SaturationModel(base_ref_files + saturation_results.meta.ref_file.saturation.name[7:])" + "saturation_ref_file = datamodels.SaturationModel(\n", + " base_ref_files + saturation_results.meta.ref_file.saturation.name[7:]\n", + ")" ] }, { @@ -993,7 +1038,7 @@ "id": "e53a8353-d769-42ec-90eb-002cf429ac24", "metadata": {}, "source": [ - "More often than not, however, the saturation reference file might not match exactly the dimensions of our subarray. This is because the reference file might be padded to match several other subarrays, and thus we have to figure out how to \"cut\" it to match our data. This is, in fact, our case:" + "More often than not, however, the saturation reference file might not match exactly the dimensions of our subarray. This is because the reference file might be padded to match several other subarrays, and thus we have to figure out how to \"cut\" it to match our data. This is, in fact, our case:\n" ] }, { @@ -1011,7 +1056,7 @@ "id": "82ebe776-8197-47d5-86ed-bcb495c328b1", "metadata": {}, "source": [ - "Luckily, the JWST calibration pipeline has a handy function to transform the dimensions between instruments --- this is the `jwst.lib.reffile_utils.get_subarray_model` function, which recieves an input data model (e.g., the one from our data) along with the reference file, and spits out the same reference file model but with the right dimensions. Let's use it:" + "Luckily, the JWST calibration pipeline has a handy function to transform the dimensions between instruments --- this is the `jwst.lib.reffile_utils.get_subarray_model` function, which recieves an input data model (e.g., the one from our data) along with the reference file, and spits out the same reference file model but with the right dimensions. Let's use it:\n" ] }, { @@ -1022,9 +1067,8 @@ "outputs": [], "source": [ "tailored_saturation_ref_file = jwst.lib.reffile_utils.get_subarray_model(\n", - " saturation_results, \n", - " saturation_ref_file\n", - " )" + " saturation_results, saturation_ref_file\n", + ")" ] }, { @@ -1032,7 +1076,7 @@ "id": "dd8b6fe1-b134-4037-a9da-16ed439ae789", "metadata": {}, "source": [ - "Indeed, now our \"tailored\" reference file matches our science data dimensions:" + "Indeed, now our \"tailored\" reference file matches our science data dimensions:\n" ] }, { @@ -1050,7 +1094,7 @@ "id": "d40b6402-c7a4-4375-b012-eb8ccce668c6", "metadata": {}, "source": [ - "Let's see how the saturation map looks like for our subarray:" + "Let's see how the saturation map looks like for our subarray:\n" ] }, { @@ -1061,9 +1105,9 @@ "outputs": [], "source": [ "plt.figure(figsize=(10, 3))\n", - "plt.title('Saturation map for NIS (SUBSSTRIP256 subarray)')\n", + "plt.title(\"Saturation map for NIS (SUBSSTRIP256 subarray)\")\n", "im = plt.imshow(tailored_saturation_ref_file.data)\n", - "plt.colorbar(label='Counts')\n", + "plt.colorbar(label=\"Counts\")\n", "plt.show()" ] }, @@ -1072,7 +1116,7 @@ "id": "98c9c2eb-cde9-49ab-b73c-d816e129131c", "metadata": {}, "source": [ - "There's clearly some structure, albeit is not exactly clear what values different pixels take. To visualize this, let's print the saturation limit for pixel `(176, 1503)`, the one we explored above:" + "There's clearly some structure, albeit is not exactly clear what values different pixels take. To visualize this, let's print the saturation limit for pixel `(176, 1503)`, the one we explored above:\n" ] }, { @@ -1083,7 +1127,7 @@ "outputs": [], "source": [ "# pixel in refernce to a saturate pixel.\n", - "tailored_saturation_ref_file.data[pixel_row, pixel_column] " + "tailored_saturation_ref_file.data[pixel_row, pixel_column]" ] }, { @@ -1091,7 +1135,7 @@ "id": "583c5962-9876-414a-8f07-69775000a6b0", "metadata": {}, "source": [ - "If the counts surpass this limit, the pixel will be considered saturated. To see if this was the case, let's repeat the plot above marking this signal limit:" + "If the counts surpass this limit, the pixel will be considered saturated. To see if this was the case, let's repeat the plot above marking this signal limit:\n" ] }, { @@ -1101,35 +1145,43 @@ "metadata": {}, "outputs": [], "source": [ - "pixel_row, pixel_column = row_idx[60], column_idx[60] \n", + "pixel_row, pixel_column = row_idx[60], column_idx[60]\n", "\n", "plt.figure(figsize=(7, 4))\n", - "plt.title(f'Saturated Pixel: ({pixel_row}, {pixel_column})')\n", - "plt.plot(np.arange(saturation_results.data.shape[1])+1, \n", - " saturation_results.data[integration, :, pixel_row, pixel_column], \n", - " 'o-', color='tomato'\n", - " )\n", + "plt.title(f\"Saturated Pixel: ({pixel_row}, {pixel_column})\")\n", + "plt.plot(\n", + " np.arange(saturation_results.data.shape[1]) + 1,\n", + " saturation_results.data[integration, :, pixel_row, pixel_column],\n", + " \"o-\",\n", + " color=\"tomato\",\n", + ")\n", "\n", - "plt.plot([1, saturation_results.data.shape[1]+1], \n", - " [tailored_saturation_ref_file.data[pixel_row, pixel_column], \n", - " tailored_saturation_ref_file.data[pixel_row, pixel_column]],\n", - " 'r--', \n", - " label='Signal limit in reference file'\n", - " )\n", + "plt.plot(\n", + " [1, saturation_results.data.shape[1] + 1],\n", + " [\n", + " tailored_saturation_ref_file.data[pixel_row, pixel_column],\n", + " tailored_saturation_ref_file.data[pixel_row, pixel_column],\n", + " ],\n", + " \"r--\",\n", + " label=\"Signal limit in reference file\",\n", + ")\n", "\n", - "plt.xlim(0.5, saturation_results.data.shape[1]+1.5)\n", - "plt.xlabel('Group number', fontsize=16)\n", - "plt.ylabel('Counts', fontsize=16, color='tomato')\n", + "plt.xlim(0.5, saturation_results.data.shape[1] + 1.5)\n", + "plt.xlabel(\"Group number\", fontsize=16)\n", + "plt.ylabel(\"Counts\", fontsize=16, color=\"tomato\")\n", "plt.legend()\n", "\n", "plt.twinx()\n", "\n", - "plt.plot(np.arange(saturation_results.data.shape[1])+1, \n", - " saturation_results.groupdq[integration, :, pixel_row, pixel_column], \n", - " 'o-', color='cornflowerblue')\n", + "plt.plot(\n", + " np.arange(saturation_results.data.shape[1]) + 1,\n", + " saturation_results.groupdq[integration, :, pixel_row, pixel_column],\n", + " \"o-\",\n", + " color=\"cornflowerblue\",\n", + ")\n", "\n", - "plt.xlim(0.5, saturation_results.data.shape[1]+1.5)\n", - "plt.ylabel('Group Data-quality', fontsize=16, color='cornflowerblue')\n", + "plt.xlim(0.5, saturation_results.data.shape[1] + 1.5)\n", + "plt.ylabel(\"Group Data-quality\", fontsize=16, color=\"cornflowerblue\")\n", "\n", "plt.show()" ] @@ -1139,7 +1191,7 @@ "id": "1a911e15-b2ae-499f-85f7-fd0e803a3834", "metadata": {}, "source": [ - "Indeed, this is the case! Note that, as described above, by default for NIRISS not only this pixel gets marked as saturated, but all pixels around it. To see this, note for instance the same plot as above but for of the neighboring pixels lets use pixel (177,1502):" + "Indeed, this is the case! Note that, as described above, by default for NIRISS not only this pixel gets marked as saturated, but all pixels around it. To see this, note for instance the same plot as above but for of the neighboring pixels lets use pixel (177,1502):\n" ] }, { @@ -1153,36 +1205,44 @@ "\n", "plt.figure(figsize=(7, 4))\n", "\n", - "plt.title(f'Same as above, but for neighboring pixel ({pixel_row},{pixel_column})')\n", - "plt.plot(np.arange(saturation_results.data.shape[1])+1, \n", - " saturation_results.data[integration, :, pixel_row, pixel_column], \n", - " 'o-', color='tomato'\n", - " )\n", - "\n", - "plt.plot([1, saturation_results.data.shape[1]+1], \n", - " [tailored_saturation_ref_file.data[pixel_row, pixel_column], \n", - " tailored_saturation_ref_file.data[pixel_row, pixel_column]],\n", - " 'r--', \n", - " label='Signal limit in reference file'\n", - " )\n", + "plt.title(f\"Same as above, but for neighboring pixel ({pixel_row},{pixel_column})\")\n", + "plt.plot(\n", + " np.arange(saturation_results.data.shape[1]) + 1,\n", + " saturation_results.data[integration, :, pixel_row, pixel_column],\n", + " \"o-\",\n", + " color=\"tomato\",\n", + ")\n", "\n", - "plt.xlim(0.5, saturation_results.data.shape[1]+1.5)\n", - "plt.xlabel('Group number', fontsize=16)\n", - "plt.ylabel('Counts', fontsize=16, color='tomato')\n", + "plt.plot(\n", + " [1, saturation_results.data.shape[1] + 1],\n", + " [\n", + " tailored_saturation_ref_file.data[pixel_row, pixel_column],\n", + " tailored_saturation_ref_file.data[pixel_row, pixel_column],\n", + " ],\n", + " \"r--\",\n", + " label=\"Signal limit in reference file\",\n", + ")\n", + "\n", + "plt.xlim(0.5, saturation_results.data.shape[1] + 1.5)\n", + "plt.xlabel(\"Group number\", fontsize=16)\n", + "plt.ylabel(\"Counts\", fontsize=16, color=\"tomato\")\n", "plt.legend()\n", "\n", "plt.twinx()\n", "\n", - "plt.plot(np.arange(saturation_results.data.shape[1])+1, \n", - " saturation_results.groupdq[integration, :, pixel_row, pixel_column], \n", - " 'o-', color='cornflowerblue')\n", + "plt.plot(\n", + " np.arange(saturation_results.data.shape[1]) + 1,\n", + " saturation_results.groupdq[integration, :, pixel_row, pixel_column],\n", + " \"o-\",\n", + " color=\"cornflowerblue\",\n", + ")\n", "\n", - "plt.xlim(0.5, saturation_results.data.shape[1]+1.5)\n", - "plt.ylabel('Group Data-quality', fontsize=16, color='cornflowerblue')\n", + "plt.xlim(0.5, saturation_results.data.shape[1] + 1.5)\n", + "plt.ylabel(\"Group Data-quality\", fontsize=16, color=\"cornflowerblue\")\n", "\n", "plt.show()\n", "\n", - "# make sure to find pixel that is saturate and its neighbor " + "# make sure to find pixel that is saturate and its neighbor" ] }, { @@ -1192,7 +1252,7 @@ "source": [ "Note how the signal level has not gone above the limit in the reference file, but it is marked as saturated because pixel (176,1503) is. Again, this is to account for possible charge spilling to the pixel.\n", "\n", - "Now, what if we wanted to mark as saturated all pixels, say, larger than 50\\% these saturation values? Well, we can directly modify the reference file and repeat the calculation pointing at it:" + "Now, what if we wanted to mark as saturated all pixels, say, larger than 50\\% these saturation values? Well, we can directly modify the reference file and repeat the calculation pointing at it:\n" ] }, { @@ -1210,7 +1270,7 @@ "id": "8a625a17-3f8d-4387-87e9-2767783dfc79", "metadata": {}, "source": [ - "To incorporate this new reference file, we simply use the `override_saturation` flag, passing this new `SaturationModel` along: " + "To incorporate this new reference file, we simply use the `override_saturation` flag, passing this new `SaturationModel` along:\n" ] }, { @@ -1222,9 +1282,8 @@ "source": [ "# Run saturation step:\n", "saturation_results2 = calwebb_detector1.saturation_step.SaturationStep.call(\n", - " uncal_nis[0], \n", - " override_saturation=saturation_ref_file\n", - " )" + " uncal_nis[0], override_saturation=saturation_ref_file\n", + ")" ] }, { @@ -1232,7 +1291,7 @@ "id": "ba32975e", "metadata": {}, "source": [ - "Let's see how many pixels are now counted as saturated:" + "Let's see how many pixels are now counted as saturated:\n" ] }, { @@ -1248,26 +1307,32 @@ "\n", "verbose = False\n", "for row in range(rows):\n", - " \n", + "\n", " for column in range(columns):\n", "\n", " # Extract the bad pixel flag(s) for the current pixel at (row, column):\n", " bps = datamodels.dqflags.dqflags_to_mnemonics(\n", - " saturation_results2.groupdq[integration, group, row, column], \n", - " mnemonic_map=datamodels.dqflags.pixel\n", + " saturation_results2.groupdq[integration, group, row, column],\n", + " mnemonic_map=datamodels.dqflags.pixel,\n", " )\n", - " \n", + "\n", " # Check if pixel is saturated; if it is...\n", - " if 'SATURATED' in bps:\n", + " if \"SATURATED\" in bps:\n", "\n", " # ...print which pixel it is, and...\n", " if verbose:\n", - " print('Pixel ({0:},{1:}) is saturated in integration 10, last group'.format(row, column))\n", + " print(\n", + " \"Pixel ({0:},{1:}) is saturated in integration 10, last group\".format(\n", + " row, column\n", + " )\n", + " )\n", "\n", " # ...count it:\n", " nsaturated += 1\n", "\n", - "print(f\"\\nA total of {nsaturated} out of {rows*columns} pixels ({100 * nsaturated / float(rows * columns):.2f}%) are saturated\")\n" + "print(\n", + " f\"\\nA total of {nsaturated} out of {rows*columns} pixels ({100 * nsaturated / float(rows * columns):.2f}%) are saturated\"\n", + ")" ] }, { @@ -1279,7 +1344,7 @@ "\n", "
Note on manually setting the saturation limit: Setting the saturation limit manually should be done with care, and we recommend trying different saturation levels to check whether TSO science is impacted by this choice. In particular, we suggest to never set limits that are above the thresholds defined by the instrument teams, as these are typically set to levels above which the non-linearity correction (see below) is not expected to work.
\n", "\n", - "Before moving to the next step, let's run the saturation step on the other NIS segments:" + "Before moving to the next step, let's run the saturation step on the other NIS segments:\n" ] }, { @@ -1302,13 +1367,13 @@ "source": [ "### 4.3 Removing detector-level effects: the `superbias` and `refpix` steps \n", "\n", - "So far, we have focused on flagging pixels for various effects (e.g., bad pixels, saturation) but we haven't worked directly with the actual counts on our data. In this Section, we deal with various (non-astrophysical) detector-level effects present in our data through two steps in the JWST Calibration pipeline: the `superbias` and the `refpix` steps. \n", + "So far, we have focused on flagging pixels for various effects (e.g., bad pixels, saturation) but we haven't worked directly with the actual counts on our data. In this Section, we deal with various (non-astrophysical) detector-level effects present in our data through two steps in the JWST Calibration pipeline: the `superbias` and the `refpix` steps.\n", "\n", "#### 4.3.1 Removing the pedestal from the detector: the `superbias` step\n", "\n", "All detectors have mostly stable, factory-defined pedestal levels, which can be closely monitored with the right calibration exposures. Indeed, instrument teams closely monitor and refine this via what is called the \"super\" bias --- the spatial shape of this pedestal. The JWST Calibration pipeline substracts this pedestal from data via the `superbias` step.\n", "\n", - "Applying this correction to the data is very simple to do; let's apply it once again to the first segment of data for the NIS detector, so we can check how our data changes after applying the step:" + "Applying this correction to the data is very simple to do; let's apply it once again to the first segment of data for the NIS detector, so we can check how our data changes after applying the step:\n" ] }, { @@ -1326,7 +1391,7 @@ "id": "6c3144d5-0b11-4984-9411-cbd4e5e87819", "metadata": {}, "source": [ - "Once again, we can see that there is a particular reference file being used to remove the pedestal, `jwst_niriss_superbias_0200.fits`, which can be explored in a similar way as how we explored the reference file for the `saturation` step above. Let's see how our data changed after applying this pedestal removal --- let's again take the last group of integration 10 as an example:" + "Once again, we can see that there is a particular reference file being used to remove the pedestal, `jwst_niriss_superbias_0200.fits`, which can be explored in a similar way as how we explored the reference file for the `saturation` step above. Let's see how our data changed after applying this pedestal removal --- let's again take the last group of integration 10 as an example:\n" ] }, { @@ -1341,19 +1406,24 @@ "ax1, ax2 = axes.ravel()\n", "\n", "# Plot before step\n", - "ax1.set_title('Before the Superbias step:')\n", - "ax1.set_ylabel('y [pixel]')\n", - "im1 = ax1.imshow(uncal_nis[0].data[10, -1, :, :] / np.nanmedian(uncal_nis[0].data[10, -1, :, :]))\n", + "ax1.set_title(\"Before the Superbias step:\")\n", + "ax1.set_ylabel(\"y [pixel]\")\n", + "im1 = ax1.imshow(\n", + " uncal_nis[0].data[10, -1, :, :] / np.nanmedian(uncal_nis[0].data[10, -1, :, :])\n", + ")\n", "im1.set_clim(-3, 2)\n", - "fig.colorbar(im1, ax=ax1, label='Normalized (to median) fluence')\n", + "fig.colorbar(im1, ax=ax1, label=\"Normalized (to median) fluence\")\n", "\n", "# Plot after step\n", - "ax2.set_title('Before the Superbias step:')\n", - "ax2.set_xlabel('x [pixel]')\n", - "ax2.set_ylabel('y [pixel]')\n", - "im2 = ax2.imshow(superbias_results.data[10, -1, :, :] / np.nanmedian(superbias_results.data[10, -1, :, :]))\n", + "ax2.set_title(\"Before the Superbias step:\")\n", + "ax2.set_xlabel(\"x [pixel]\")\n", + "ax2.set_ylabel(\"y [pixel]\")\n", + "im2 = ax2.imshow(\n", + " superbias_results.data[10, -1, :, :]\n", + " / np.nanmedian(superbias_results.data[10, -1, :, :])\n", + ")\n", "im2.set_clim(-3, 2)\n", - "fig.colorbar(im2, ax=ax2, label='Normalized (to median) fluence')\n", + "fig.colorbar(im2, ax=ax2, label=\"Normalized (to median) fluence\")\n", "\n", "plt.tight_layout()\n", "plt.show()" @@ -1364,7 +1434,7 @@ "id": "340e93a6", "metadata": {}, "source": [ - "Wow! That's a huge change. Overall, this looks much better and the 3rd diffraction order that was previously being suppressed by detector effect is now visible. Let's plot the profiles of pixel column index 1500 to have a closer look:" + "Wow! That's a huge change. Overall, this looks much better and the 3rd diffraction order that was previously being suppressed by detector effect is now visible. Let's plot the profiles of pixel column index 1500 to have a closer look:\n" ] }, { @@ -1380,18 +1450,26 @@ "\n", "ax1, ax2 = axes.ravel()\n", "\n", - "ax1.plot(uncal_nis[0].data[10, -1, :, column_index], label='Before the Superbias step')\n", - "ax1.plot(superbias_results.data[10, -1, :, column_index], label='After the Superbias step')\n", - "ax1.set_xlabel('Row pixel index', fontsize=14)\n", - "ax1.set_ylabel('Counts', fontsize=14)\n", - "ax1.set_title('Comparison before/after Superbias step', fontsize=14)\n", + "ax1.plot(uncal_nis[0].data[10, -1, :, column_index], label=\"Before the Superbias step\")\n", + "ax1.plot(\n", + " superbias_results.data[10, -1, :, column_index], label=\"After the Superbias step\"\n", + ")\n", + "ax1.set_xlabel(\"Row pixel index\", fontsize=14)\n", + "ax1.set_ylabel(\"Counts\", fontsize=14)\n", + "ax1.set_title(\"Comparison before/after Superbias step\", fontsize=14)\n", "ax1.legend()\n", "\n", - "ax2.set_title('Same, but median-subtracted counts', fontsize=14)\n", - "ax2.plot(uncal_nis[0].data[10, -1, :, column_index] - np.nanmedian(uncal_nis[0].data[0, -1, :, column_index]))\n", - "ax2.plot(superbias_results.data[10, -1, :, column_index] - np.nanmedian(superbias_results.data[0, -1, :, column_index]))\n", - "ax2.set_xlabel('Row pixel index', fontsize=14)\n", - "ax2.set_ylabel('Counts - Median Counts', fontsize=14)\n", + "ax2.set_title(\"Same, but median-subtracted counts\", fontsize=14)\n", + "ax2.plot(\n", + " uncal_nis[0].data[10, -1, :, column_index]\n", + " - np.nanmedian(uncal_nis[0].data[0, -1, :, column_index])\n", + ")\n", + "ax2.plot(\n", + " superbias_results.data[10, -1, :, column_index]\n", + " - np.nanmedian(superbias_results.data[0, -1, :, column_index])\n", + ")\n", + "ax2.set_xlabel(\"Row pixel index\", fontsize=14)\n", + "ax2.set_ylabel(\"Counts - Median Counts\", fontsize=14)\n", "\n", "plt.show()" ] @@ -1401,9 +1479,9 @@ "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, in principle, correctly. \n", + "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, in principle, correctly.\n", "\n", - "However, if we take a more detailed look at background pixels, we can note an interesting pattern. Let's plot a similar cut to the one above, but for column 250 --- which is far away from any illuminted pixels in the detector. Let's also plot the last superbias-corrected group and the second-to-last superbias-corrected group:" + "However, if we take a more detailed look at background pixels, we can note an interesting pattern. Let's plot a similar cut to the one above, but for column 250 --- which is far away from any illuminted pixels in the detector. Let's also plot the last superbias-corrected group and the second-to-last superbias-corrected group:\n" ] }, { @@ -1418,20 +1496,22 @@ "ax1, ax2 = axes.ravel()\n", "\n", "ax1.plot(superbias_results.data[10, -1, :, 250])\n", - "ax1.plot([0, 32], [0, 0], 'k--')\n", - "ax1.set_xlabel('Row pixel index', fontsize=14)\n", - "ax1.set_ylabel('Counts', fontsize=14)\n", - "ax1.set_title('Superbias-corrected close-up, last group, integration 10', fontsize=14)\n", + "ax1.plot([0, 32], [0, 0], \"k--\")\n", + "ax1.set_xlabel(\"Row pixel index\", fontsize=14)\n", + "ax1.set_ylabel(\"Counts\", fontsize=14)\n", + "ax1.set_title(\"Superbias-corrected close-up, last group, integration 10\", fontsize=14)\n", "ax1.set_ylim(-250, 250)\n", "ax1.set_xlim(0, 31)\n", "\n", - "ax2.set_title('Superbias-corrected close-up, second-to-last group, integration 10', fontsize=14)\n", + "ax2.set_title(\n", + " \"Superbias-corrected close-up, second-to-last group, integration 10\", fontsize=14\n", + ")\n", "ax2.plot(superbias_results.data[10, -2, :, 250])\n", - "ax2.plot([0, 32], [0, 0], 'k--')\n", + "ax2.plot([0, 32], [0, 0], \"k--\")\n", "ax2.set_ylim(-250, 250)\n", "ax2.set_xlim(0, 31)\n", - "ax2.set_xlabel('Row pixel index', fontsize=14)\n", - "ax2.set_ylabel('Counts', fontsize=14)\n", + "ax2.set_xlabel(\"Row pixel index\", fontsize=14)\n", + "ax2.set_ylabel(\"Counts\", fontsize=14)\n", "\n", "plt.show()" ] @@ -1447,7 +1527,7 @@ "\n", "All the JWST detectors contain reference pixels, typically located in some (or all) of the edges of the detectors. These pixels are ones for which their \"sensitivity to light\" has been deactivated, and are thus useful for tracking detector-level effects happening at the time of our observations. While all detectors have those, **not all detector subarrays** contain reference pixels. Some, like in our case, contain reference pixels only in certain portions of the subarray.\n", "\n", - "Let's visualize where those reference pixels are in our subarray by using the `pixeldq` flags:" + "Let's visualize where those reference pixels are in our subarray by using the `pixeldq` flags:\n" ] }, { @@ -1462,16 +1542,16 @@ "\n", "# Iterate through every row and column:\n", "for row in range(rows):\n", - " \n", + "\n", " for column in range(columns):\n", "\n", " # Extract the bad pixel flag(s) for the current pixel at (row, column):\n", " bps = datamodels.dqflags.dqflags_to_mnemonics(\n", - " superbias_results.pixeldq[row, column], \n", - " mnemonic_map=datamodels.dqflags.pixel\n", + " superbias_results.pixeldq[row, column],\n", + " mnemonic_map=datamodels.dqflags.pixel,\n", " )\n", "\n", - " if 'REFERENCE_PIXEL' in bps:\n", + " if \"REFERENCE_PIXEL\" in bps:\n", "\n", " reference_pixels[row, column] = 1" ] @@ -1485,19 +1565,28 @@ "source": [ "plt.figure(figsize=(12, 4))\n", "\n", - "plt.title('Location of reference pixels in the subarray: 4 columns/rows')\n", + "plt.title(\"Location of reference pixels in the subarray: 4 columns/rows\")\n", "im = plt.imshow(reference_pixels)\n", "\n", "# Arrows to indicate edges:\n", - "plt.text(1800-70, 32, 'Right Edge Ref Pixels', color='yellow')\n", - "plt.arrow(1780, 16, 150, -1, widt=5, head_width=10, head_length=100, color='white')\n", + "plt.text(1800 - 70, 32, \"Right Edge Ref Pixels\", color=\"yellow\")\n", + "plt.arrow(1780, 16, 150, -1, widt=5, head_width=10, head_length=100, color=\"white\")\n", "\n", - "plt.text(30, 32, 'Left Edge Ref Pixels', color='yellow')\n", - "plt.arrow(268, 16, -150, -1, width=5, head_width=10, head_length=100, color='white')\n", + "plt.text(30, 32, \"Left Edge Ref Pixels\", color=\"yellow\")\n", + "plt.arrow(268, 16, -150, -1, width=5, head_width=10, head_length=100, color=\"white\")\n", "\n", "# plot vertical error top reference pixels\n", - "plt.text(1044, 170, 'Top Ref Pixels', color='yellow', rotation=90)\n", - "plt.arrow(1024, 180, -0, 150//4, width=5*2, head_width=10*2, head_length=100/2/2, color='white')\n", + "plt.text(1044, 170, \"Top Ref Pixels\", color=\"yellow\", rotation=90)\n", + "plt.arrow(\n", + " 1024,\n", + " 180,\n", + " -0,\n", + " 150 // 4,\n", + " width=5 * 2,\n", + " head_width=10 * 2,\n", + " head_length=100 / 2 / 2,\n", + " color=\"white\",\n", + ")\n", "\n", "im.set_clim(-0.5, 0.5)\n", "\n", @@ -1509,7 +1598,7 @@ "id": "a183faf7-693b-4ca0-97da-321d9703012a", "metadata": {}, "source": [ - "Note the white/yellow edges to the left and the right of the plot above, indicated by arrows, show the location of the reference pixels for our subarray. In other words, our subarray has reference pixels at the top of the frame andto the left and right-most sides, but not the bottom part." + "Note the white/yellow edges to the left and the right of the plot above, indicated by arrows, show the location of the reference pixels for our subarray. In other words, our subarray has reference pixels at the top of the frame andto the left and right-most sides, but not the bottom part.\n" ] }, { @@ -1517,7 +1606,7 @@ "id": "67c3e82f", "metadata": {}, "source": [ - "Let's apply the `refpix` step to check how our data looks like after it:" + "Let's apply the `refpix` step to check how our data looks like after it:\n" ] }, { @@ -1537,7 +1626,7 @@ "source": [ "Let's plot the before and after applying this step:\n", "\n", - "" + "\n" ] }, { @@ -1552,26 +1641,31 @@ "ax1, ax2, ax3 = axes.ravel()\n", "\n", "# Plot before step\n", - "ax1.set_title('Before the RefPix step:')\n", - "ax1.set_ylabel('y [pixel]')\n", - "im1 = ax1.imshow(superbias_results.data[0, -1, :, :] / np.nanmedian(superbias_results.data[0, -1, :, :]))\n", + "ax1.set_title(\"Before the RefPix step:\")\n", + "ax1.set_ylabel(\"y [pixel]\")\n", + "im1 = ax1.imshow(\n", + " superbias_results.data[0, -1, :, :]\n", + " / np.nanmedian(superbias_results.data[0, -1, :, :])\n", + ")\n", "im1.set_clim(-2, 2)\n", - "fig.colorbar(im1, ax=ax1, label='Normalized (to median) fluence')\n", + "fig.colorbar(im1, ax=ax1, label=\"Normalized (to median) fluence\")\n", "\n", "# Plot after step\n", - "ax2.set_title('Before the RefPix step:')\n", - "ax2.set_xlabel('x [pixel]')\n", - "ax2.set_ylabel('y [pixel]')\n", - "im2 = ax2.imshow(refpix_results.data[0, -1, :, :] / np.nanmedian(refpix_results.data[0, -1, :, :]))\n", + "ax2.set_title(\"Before the RefPix step:\")\n", + "ax2.set_xlabel(\"x [pixel]\")\n", + "ax2.set_ylabel(\"y [pixel]\")\n", + "im2 = ax2.imshow(\n", + " refpix_results.data[0, -1, :, :] / np.nanmedian(refpix_results.data[0, -1, :, :])\n", + ")\n", "im2.set_clim(-2, 2)\n", - "fig.colorbar(im2, ax=ax2, label='Normalized (to median) fluence')\n", + "fig.colorbar(im2, ax=ax2, label=\"Normalized (to median) fluence\")\n", "\n", "# Plot difference\n", - "ax3.set_title('Difference:')\n", - "ax3.set_xlabel('x [pixel]')\n", - "ax3.set_ylabel('y [pixel]')\n", + "ax3.set_title(\"Difference:\")\n", + "ax3.set_xlabel(\"x [pixel]\")\n", + "ax3.set_ylabel(\"y [pixel]\")\n", "im3 = ax3.imshow(superbias_results.data[0, -1, :, :] - refpix_results.data[0, -1, :, :])\n", - "fig.colorbar(im3, ax=ax3, label='Normalized (to median) fluence')\n", + "fig.colorbar(im3, ax=ax3, label=\"Normalized (to median) fluence\")\n", "\n", "plt.tight_layout()\n", "plt.show()" @@ -1584,7 +1678,7 @@ "source": [ "That looks much better. Several things have been removed, including some interesting high-frequency noise happening in the rows as show in the difference between the before and after. That's the so-called odd-even effect, 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." + "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" ] }, { @@ -1598,19 +1692,25 @@ "\n", "ax1, ax2 = axes.ravel()\n", "\n", - "ax1.plot(superbias_results.data[0, -1, :, column_index], label='Before the RefPix step')\n", - "ax1.plot(refpix_results.data[0, -1, :, column_index], label='After the RefPix step')\n", - "ax1.set_xlabel('Row pixel index', fontsize=14)\n", - "ax1.set_ylabel('Counts', fontsize=14)\n", - "ax1.set_title('Comparison before/after RefPix step', fontsize=14)\n", + "ax1.plot(superbias_results.data[0, -1, :, column_index], label=\"Before the RefPix step\")\n", + "ax1.plot(refpix_results.data[0, -1, :, column_index], label=\"After the RefPix step\")\n", + "ax1.set_xlabel(\"Row pixel index\", fontsize=14)\n", + "ax1.set_ylabel(\"Counts\", fontsize=14)\n", + "ax1.set_title(\"Comparison before/after RefPix step\", fontsize=14)\n", "ax1.legend()\n", "\n", "\n", - "ax2.set_title('Same, but median-substracted counts', fontsize=14)\n", - "ax2.plot(superbias_results.data[0, -1, :, column_index] - np.nanmedian(superbias_results.data[0, -1, :, column_index]))\n", - "ax2.plot(refpix_results.data[0, -1, :, column_index] - np.nanmedian(refpix_results.data[0, -1, :, column_index]))\n", - "ax2.set_xlabel('Row pixel index', fontsize=14)\n", - "ax2.set_ylabel('Counts - Median Counts', fontsize=14)\n", + "ax2.set_title(\"Same, but median-substracted counts\", fontsize=14)\n", + "ax2.plot(\n", + " superbias_results.data[0, -1, :, column_index]\n", + " - np.nanmedian(superbias_results.data[0, -1, :, column_index])\n", + ")\n", + "ax2.plot(\n", + " refpix_results.data[0, -1, :, column_index]\n", + " - np.nanmedian(refpix_results.data[0, -1, :, column_index])\n", + ")\n", + "ax2.set_xlabel(\"Row pixel index\", fontsize=14)\n", + "ax2.set_ylabel(\"Counts - Median Counts\", fontsize=14)\n", "\n", "plt.tight_layout()\n", "plt.show()" @@ -1621,9 +1721,9 @@ "id": "7cdf297f-1f61-4c05-a6c6-d00909136f3c", "metadata": {}, "source": [ - "We can see there's been a minor improvement in the results. \n", + "We can see there's been a minor improvement in the results.\n", "\n", - "It is interesting to note that the \"banding\" on the columns, as discussed above, has not dissapeared. This is more evident when plotting a series of groups from different integrations; let's plots the groups from integrations 10, 11 and 12:" + "It is interesting to note that the \"banding\" on the columns, as discussed above, has not dissapeared. This is more evident when plotting a series of groups from different integrations; let's plots the groups from integrations 10, 11 and 12:\n" ] }, { @@ -1635,24 +1735,24 @@ "source": [ "# Integration 10, last group\n", "plt.figure(figsize=(10, 3))\n", - "plt.title('Integration 10, last group')\n", + "plt.title(\"Integration 10, last group\")\n", "im = plt.imshow(refpix_results.data[10, -1, :, :])\n", "im.set_clim(-200, 400)\n", - "plt.colorbar(label='Counts')\n", + "plt.colorbar(label=\"Counts\")\n", "\n", "# Integration 10, last group\n", "plt.figure(figsize=(10, 3))\n", - "plt.title('Integration 11, last group')\n", + "plt.title(\"Integration 11, last group\")\n", "im = plt.imshow(refpix_results.data[11, -1, :, :])\n", "im.set_clim(-200, 400)\n", - "plt.colorbar(label='Counts')\n", + "plt.colorbar(label=\"Counts\")\n", "\n", "# Integration 12, last group\n", "plt.figure(figsize=(10, 3))\n", - "plt.title('Integration 12, last group')\n", + "plt.title(\"Integration 12, last group\")\n", "im = plt.imshow(refpix_results.data[12, -1, :, :])\n", "im.set_clim(-200, 400)\n", - "plt.colorbar(label='Counts')\n", + "plt.colorbar(label=\"Counts\")\n", "\n", "plt.show()" ] @@ -1662,7 +1762,7 @@ "id": "3837a984-6b6e-4583-9077-1a5456a480c3", "metadata": {}, "source": [ - "This is, once again, expected as there are no reference pixels in the columns. We will explore how to correct this after going with the `linearity` correction/step, which we discuss next. Before moving on, we apply the superbias and reference pixel step to both detectors, all segments:" + "This is, once again, expected as there are no reference pixels in the columns. We will explore how to correct this after going with the `linearity` correction/step, which we discuss next. Before moving on, we apply the superbias and reference pixel step to both detectors, all segments:\n" ] }, { @@ -1691,7 +1791,7 @@ "\n", "#### 4.4.1 Visualizing and correcting for non-linearity with the `linearity` step\n", "\n", - "To visualize the non-linearity of the up-the-ramp samples, let's take a look at the samples of one of the brightest pixels in our subarray, pixel `(45, 1600)` --- say for integration number 10. Let's plot on top a line fitted to the first 10 pixels, which should be the most \"linear\" of all pixels:" + "To visualize the non-linearity of the up-the-ramp samples, let's take a look at the samples of one of the brightest pixels in our subarray, pixel `(45, 1600)` --- say for integration number 10. Let's plot on top a line fitted to the first 10 pixels, which should be the most \"linear\" of all pixels:\n" ] }, { @@ -1707,24 +1807,38 @@ "first_groups = 9\n", "i1, i2 = 45, 1600\n", "\n", - "coeff = np.polyfit(group[:first_groups], uncal_nis[0].data[10, :first_groups, i1, i2], 1)\n", + "coeff = np.polyfit(\n", + " group[:first_groups], uncal_nis[0].data[10, :first_groups, i1, i2], 1\n", + ")\n", "\n", "fig, axes = plt.subplots(2, 1, figsize=(6, 6), constrained_layout=True)\n", "\n", "ax1, ax2 = axes.ravel()\n", "\n", "# plot ramp samples\n", - "ax1.set_title(f'Up-the-ramp sample, integration 10, pixel ({i1}, {i2})')\n", - "ax1.plot(group, uncal_nis[0].data[10, :, i1, i2], 'o-', color='black', mfc='white', label='Up-the-ramp samples')\n", - "ax1.plot(group, np.polyval(coeff, group), 'r--', label='Linear fit to first ' + str(first_groups) + ' groups')\n", - "ax1.set_xlabel('Group number', fontsize=16)\n", - "ax1.set_ylabel('Counts', fontsize=16)\n", + "ax1.set_title(f\"Up-the-ramp sample, integration 10, pixel ({i1}, {i2})\")\n", + "ax1.plot(\n", + " group,\n", + " uncal_nis[0].data[10, :, i1, i2],\n", + " \"o-\",\n", + " color=\"black\",\n", + " mfc=\"white\",\n", + " label=\"Up-the-ramp samples\",\n", + ")\n", + "ax1.plot(\n", + " group,\n", + " np.polyval(coeff, group),\n", + " \"r--\",\n", + " label=\"Linear fit to first \" + str(first_groups) + \" groups\",\n", + ")\n", + "ax1.set_xlabel(\"Group number\", fontsize=16)\n", + "ax1.set_ylabel(\"Counts\", fontsize=16)\n", "ax1.legend()\n", "ax1.set_xlim(0.5, 9.5)\n", "\n", "# plot residuals\n", - "ax2.plot(group, uncal_nis[0].data[10, :, i1, i2] - np.polyval(coeff, group), 'o-')\n", - "ax2.set_xlabel('Group number', fontsize=16)\n", + "ax2.plot(group, uncal_nis[0].data[10, :, i1, i2] - np.polyval(coeff, group), \"o-\")\n", + "ax2.set_xlabel(\"Group number\", fontsize=16)\n", "ax2.set_ylabel(\"residuals\")\n", "plt.show()" ] @@ -1734,7 +1848,7 @@ "id": "7d261731-04de-466c-8f0e-e9dea55af661", "metadata": {}, "source": [ - "Ah --- the ramp is _clearly_ non-linear! Let's apply the `linearity` step to the very first segment to see how well this gets corrected:" + "Ah --- the ramp is _clearly_ non-linear! Let's apply the `linearity` step to the very first segment to see how well this gets corrected:\n" ] }, { @@ -1753,7 +1867,7 @@ "id": "2882aaf3-0f2d-492b-9a97-07e06ca44d72", "metadata": {}, "source": [ - "Let's try the same plot as above, but with the linearity-corrected data:" + "Let's try the same plot as above, but with the linearity-corrected data:\n" ] }, { @@ -1768,24 +1882,38 @@ "\n", "first_groups = 9\n", "\n", - "coeff = np.polyfit(group[:first_groups], linearity_results.data[10, :first_groups, i1, i2], 1)\n", + "coeff = np.polyfit(\n", + " group[:first_groups], linearity_results.data[10, :first_groups, i1, i2], 1\n", + ")\n", "\n", "fig, axes = plt.subplots(2, 1, figsize=(6, 6), constrained_layout=True)\n", "\n", "ax1, ax2 = axes.ravel()\n", "\n", "# plot ramp samples\n", - "ax1.set_title('Same as above, linearity-corrected')\n", - "ax1.plot(group, linearity_results.data[10, :, i1, i2], 'o-', color='black', mfc='white', label='Up-the-ramp samples')\n", - "ax1.plot(group, np.polyval(coeff, group), 'r--', label='Linear fit to first ' + str(first_groups) + ' groups')\n", - "ax1.set_xlabel('Group number', fontsize=16)\n", - "ax1.set_ylabel('Counts', fontsize=16)\n", + "ax1.set_title(\"Same as above, linearity-corrected\")\n", + "ax1.plot(\n", + " group,\n", + " linearity_results.data[10, :, i1, i2],\n", + " \"o-\",\n", + " color=\"black\",\n", + " mfc=\"white\",\n", + " label=\"Up-the-ramp samples\",\n", + ")\n", + "ax1.plot(\n", + " group,\n", + " np.polyval(coeff, group),\n", + " \"r--\",\n", + " label=\"Linear fit to first \" + str(first_groups) + \" groups\",\n", + ")\n", + "ax1.set_xlabel(\"Group number\", fontsize=16)\n", + "ax1.set_ylabel(\"Counts\", fontsize=16)\n", "ax1.legend()\n", "ax1.set_xlim(0.5, 9.5)\n", "\n", "# plot residuals\n", - "ax2.plot(group, linearity_results.data[10, :, i1, i2] - np.polyval(coeff, group), 'o-')\n", - "ax2.set_xlabel('Group number', fontsize=16)\n", + "ax2.plot(group, linearity_results.data[10, :, i1, i2] - np.polyval(coeff, group), \"o-\")\n", + "ax2.set_xlabel(\"Group number\", fontsize=16)\n", "ax2.set_ylabel(\"residuals\")\n", "plt.show()" ] @@ -1795,11 +1923,11 @@ "id": "6ca8a4da-90ef-4581-a06e-f8a8395e334a", "metadata": {}, "source": [ - "Ah, much better, notice how the residuals are have improve and are much closer to zero. \n", + "Ah, much better, notice how the residuals are have improve and are much closer to zero.\n", "\n", "#### 4.4.2 Testing the accuracy of the `linearity` step\n", "\n", - "It is important to realize that the linearity corrections that the JWST Calibration pipeline applies through the `linearity` step are _not_ perfect. While this is difficult to see with a single integration, this can be studied with multiple integrations --- which helps us beat the noise embedded on single up-the-ramp samples. \n", + "It is important to realize that the linearity corrections that the JWST Calibration pipeline applies through the `linearity` step are _not_ perfect. While this is difficult to see with a single integration, this can be studied with multiple integrations --- which helps us beat the noise embedded on single up-the-ramp samples.\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", @@ -1807,7 +1935,7 @@ "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 of the **uncorrected** data --- then take the ratio of those. As we saw above, this should scream \"non-linearity\" all over!" + "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 of the **uncorrected** data --- then take the ratio of those. As we saw above, this should scream \"non-linearity\" all over!\n" ] }, { @@ -1836,7 +1964,7 @@ "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 (this may take a few seconds):" + "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 (this may take a few seconds):\n" ] }, { @@ -1858,13 +1986,13 @@ "outputs": [], "source": [ "plt.figure(figsize=(6, 4))\n", - "plt.plot(flattened_fluences, flattened_ratio, '.', alpha=0.01, color='black')\n", - "plt.plot([0, 35000], [1., 1.], 'r--')\n", + "plt.plot(flattened_fluences, flattened_ratio, \".\", alpha=0.01, color=\"black\")\n", + "plt.plot([0, 35000], [1.0, 1.0], \"r--\")\n", "plt.ylim(0.5, 1.5)\n", "plt.xlim(0, 35000)\n", - "plt.xlabel('Fluence at the last group (counts)', fontsize=14)\n", - "plt.ylabel('(Last / First) Group differences', fontsize=14)\n", - "plt.title('No linearity correction', fontsize=14)\n", + "plt.xlabel(\"Fluence at the last group (counts)\", fontsize=14)\n", + "plt.ylabel(\"(Last / First) Group differences\", fontsize=14)\n", + "plt.title(\"No linearity correction\", fontsize=14)\n", "\n", "plt.show()" ] @@ -1876,7 +2004,7 @@ "source": [ "Indeed, the data does _not_ line up around 1. So linearity _is_ an issue the larger the flux received (as we already observed in the up-the ramp samples before)!\n", "\n", - "All right, let's try the same experiment but now on the linearity-corrected data:" + "All right, let's try the same experiment but now on the linearity-corrected data:\n" ] }, { @@ -1886,8 +2014,12 @@ "metadata": {}, "outputs": [], "source": [ - "corrected_last_pair = linearity_results.data[:, -1, :, :] - linearity_results.data[:, -2, :, :]\n", - "corrected_first_pair = linearity_results.data[:, 1, :, :] - linearity_results.data[:, 0, :, :]\n", + "corrected_last_pair = (\n", + " linearity_results.data[:, -1, :, :] - linearity_results.data[:, -2, :, :]\n", + ")\n", + "corrected_first_pair = (\n", + " linearity_results.data[:, 1, :, :] - linearity_results.data[:, 0, :, :]\n", + ")\n", "corrected_ratio = corrected_last_pair / corrected_first_pair" ] }, @@ -1896,7 +2028,7 @@ "id": "3041b82e", "metadata": {}, "source": [ - "Let's plot:" + "Let's plot:\n" ] }, { @@ -1918,13 +2050,19 @@ "outputs": [], "source": [ "plt.figure(figsize=(6, 4))\n", - "plt.plot(flattened_corrected_fluences, flattened_corrected_ratio, '.', alpha=0.005, color='black')\n", - "plt.plot([0, 35000], [1., 1.], 'r--')\n", + "plt.plot(\n", + " flattened_corrected_fluences,\n", + " flattened_corrected_ratio,\n", + " \".\",\n", + " alpha=0.005,\n", + " color=\"black\",\n", + ")\n", + "plt.plot([0, 35000], [1.0, 1.0], \"r--\")\n", "plt.ylim(0.5, 1.5)\n", "plt.xlim(0, 35000)\n", - "plt.xlabel('Fluence at the last group (counts)', fontsize=14)\n", - "plt.ylabel('(Last / First) Group differences', fontsize=14)\n", - "plt.title('After linearity correction', fontsize=14)\n", + "plt.xlabel(\"Fluence at the last group (counts)\", fontsize=14)\n", + "plt.ylabel(\"(Last / First) Group differences\", fontsize=14)\n", + "plt.title(\"After linearity correction\", fontsize=14)\n", "plt.show()" ] }, @@ -1933,9 +2071,9 @@ "id": "1896c9fd-da47-4def-9882-62cde592a8bb", "metadata": {}, "source": [ - "That looks **much** better. Note, however, that as discussed above the corrections are *not* perfect. In particular, below about 20,000 counts it seems the correction makes the last group difference to be slightly larger than the first group differences; this changes for the larger fluences, where the last group difference seems to have a _smaller_ flux than the first group differences. This is actually consistent with a _charge migration_ hypothesis, on which pixels that receive larger fluences _lose_ charge to neighboring pixels that receive them. Testing this hypothesis is, of course, outside of the present Notebook --- but this showcases that plots like the ones above are fundamental to make sense of data and the overall accuracy and precision of non-linearity corrections.\n", + "That looks **much** better. Note, however, that as discussed above the corrections are _not_ perfect. In particular, below about 20,000 counts it seems the correction makes the last group difference to be slightly larger than the first group differences; this changes for the larger fluences, where the last group difference seems to have a _smaller_ flux than the first group differences. This is actually consistent with a _charge migration_ hypothesis, on which pixels that receive larger fluences _lose_ charge to neighboring pixels that receive them. Testing this hypothesis is, of course, outside of the present Notebook --- but this showcases that plots like the ones above are fundamental to make sense of data and the overall accuracy and precision of non-linearity corrections.\n", "\n", - "Before moving to the next step, we apply the `linearity` step to all our data:" + "Before moving to the next step, we apply the `linearity` step to all our data:\n" ] }, { @@ -1959,9 +2097,9 @@ "source": [ "### 4.5 Removing the Dark Current \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. \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.\n", "\n", - "It is unclear if this step is helpful at all for TSOs, where signals are typically high (and thus, the dark current is but a very small addition to the total current gathered in a TSO), but we go ahead and apply this step nonetheless in our data. First, to check what changes this step does in our data, we apply it on the first NIS segment: " + "It is unclear if this step is helpful at all for TSOs, where signals are typically high (and thus, the dark current is but a very small addition to the total current gathered in a TSO), but we go ahead and apply this step nonetheless in our data. First, to check what changes this step does in our data, we apply it on the first NIS segment:\n" ] }, { @@ -1972,7 +2110,9 @@ "outputs": [], "source": [ "# Run the darkcurrent step:\n", - "darkcurrent_results = calwebb_detector1.dark_current_step.DarkCurrentStep.call(uncal_nis[0])" + "darkcurrent_results = calwebb_detector1.dark_current_step.DarkCurrentStep.call(\n", + " uncal_nis[0]\n", + ")" ] }, { @@ -1980,7 +2120,7 @@ "id": "755efd9f", "metadata": {}, "source": [ - "Let's see its impact on products before the dark current correction:" + "Let's see its impact on products before the dark current correction:\n" ] }, { @@ -1992,17 +2132,22 @@ "source": [ "# Plot them:\n", "plt.figure(figsize=(10, 3))\n", - "im = plt.imshow(uncal_nis[0].data[10, -1, :, :] / np.nanmedian(uncal_nis[0].data[10, -1, :, :]))\n", + "im = plt.imshow(\n", + " uncal_nis[0].data[10, -1, :, :] / np.nanmedian(uncal_nis[0].data[10, -1, :, :])\n", + ")\n", "im.set_clim(-3, 2)\n", - "plt.colorbar(label='Normalized (to median) fluence')\n", - "plt.title('Before the DarkCurrent step:')\n", + "plt.colorbar(label=\"Normalized (to median) fluence\")\n", + "plt.title(\"Before the DarkCurrent step:\")\n", "\n", "# Plot them:\n", "plt.figure(figsize=(10, 3))\n", - "im = plt.imshow(darkcurrent_results.data[10, -1, :, :] / np.nanmedian(darkcurrent_results.data[10, -1, :, :]))\n", + "im = plt.imshow(\n", + " darkcurrent_results.data[10, -1, :, :]\n", + " / np.nanmedian(darkcurrent_results.data[10, -1, :, :])\n", + ")\n", "im.set_clim(-3, 2)\n", - "plt.title('After the DarkCurrent step:')\n", - "plt.colorbar(label='Normalized (to median) fluence')\n", + "plt.title(\"After the DarkCurrent step:\")\n", + "plt.colorbar(label=\"Normalized (to median) fluence\")\n", "\n", "plt.show()" ] @@ -2012,11 +2157,11 @@ "id": "342f349d", "metadata": {}, "source": [ - "Difficult to see the impact from those simple plots. \n", + "Difficult to see the impact from those simple plots.\n", "\n", - "Let's quantify \"how much\" dark current there is by simply calculating the average (accross integrations) percentage of dark current on the last group. The reason for doing this in the last group is that this is the group that accumulates _the most_ dark current --- as dark current grows as a function of the number of groups. \n", + "Let's quantify \"how much\" dark current there is by simply calculating the average (accross integrations) percentage of dark current on the last group. The reason for doing this in the last group is that this is the group that accumulates _the most_ dark current --- as dark current grows as a function of the number of groups.\n", "\n", - "To do this, we consider that for a non-dark current corrected signal $S_{DC}$, if we substract the dark-current corrected signal $S_{DC, corrected}$ we get the dark current signal back, i.e., $S_{DC} - S_{DC, corrected} = DC$; dividing this by the dark-current corrected signal gives us the percentage of dark current signal on each pixel. Let's calculate a map of this for all integrations and take the median of those:" + "To do this, we consider that for a non-dark current corrected signal $S_{DC}$, if we substract the dark-current corrected signal $S_{DC, corrected}$ we get the dark current signal back, i.e., $S_{DC} - S_{DC, corrected} = DC$; dividing this by the dark-current corrected signal gives us the percentage of dark current signal on each pixel. Let's calculate a map of this for all integrations and take the median of those:\n" ] }, { @@ -2027,8 +2172,10 @@ "outputs": [], "source": [ "# Calculate the (average) percent change of the signal accross all integrations --- this is (Dark Signal) / (\"Real\" signal):\n", - "percent = ((linearity_results.data[:, -1, :, :] - darkcurrent_results.data[:, -1, :, :]) / \n", - " darkcurrent_results.data[:, -1, :, :]) * 100\n", + "percent = (\n", + " (linearity_results.data[:, -1, :, :] - darkcurrent_results.data[:, -1, :, :])\n", + " / darkcurrent_results.data[:, -1, :, :]\n", + ") * 100\n", "\n", "percent = np.nanmedian(percent, axis=0)\n", "\n", @@ -2036,8 +2183,8 @@ "plt.figure(figsize=(10, 3))\n", "im = plt.imshow(percent)\n", "im.set_clim(0, 25)\n", - "plt.colorbar(label='% of Dark Signal')\n", - "plt.title('Median impact of dark signal on the last group')\n", + "plt.colorbar(label=\"% of Dark Signal\")\n", + "plt.title(\"Median impact of dark signal on the last group\")\n", "plt.show()" ] }, @@ -2050,7 +2197,7 @@ "\n", "How much the above impacts a given TSO must be defined on a target-by-target basis. In the worst-case scenarios, this impact might not be purely aesthetical --- this dark current can give rise to transit depth dilutions in transiting exoplanet science, for instance; just as any non-accounted background signal.\n", "\n", - "In the case of this notebook, we apply it nonetheless to all the detector-level data in all segments, but we leave as an excercise to the reader to perform a full re-reduction with and without dark-current to see the impact of this step on this particular dataset:" + "In the case of this notebook, we apply it nonetheless to all the detector-level data in all segments, but we leave as an excercise to the reader to perform a full re-reduction with and without dark-current to see the impact of this step on this particular dataset:\n" ] }, { @@ -2064,7 +2211,9 @@ "source": [ "for i in range(nsegments):\n", " # Apply the dark_current step to NIS segments:\n", - " uncal_nis[i] = calwebb_detector1.dark_current_step.DarkCurrentStep.call(uncal_nis[i])" + " uncal_nis[i] = calwebb_detector1.dark_current_step.DarkCurrentStep.call(\n", + " uncal_nis[i]\n", + " )" ] }, { @@ -2076,7 +2225,7 @@ "\n", "For TSOs, there is some discussion in the literature about whether attempting to remove 1/f noise at the group-level, at the rate-level (i.e., after all the steps in detector1) or both is the way to go --- and whether simplistic algorithms provide a quick means of removed this source of noise. The reality is that, at the time of writing, the jury is still out on the final answer. We thus encourage readers to try different methodologies and find the one that works best for their scientific use-case. As a start, an interesting reader might, e.g., skip the above 1/f removal algorithm and simply try to remove it at the rate-level --- or perform no removal at all, and see differences in the final lightcurve precision.\n", "\n", - "" + "\n" ] }, { @@ -2086,13 +2235,12 @@ "source": [ "### 4.7 Detecting \"jumps\" in up-the-ramp samples \n", "\n", - "When a cosmic-ray hits JWST detectors, this impacts the up-the-ramp samples by making them \"[jump](https://www.youtube.com/watch?v=SwYN7mTi6HM)\" from one group to another. We already noted this happening above \n", + "When a cosmic-ray hits JWST detectors, this impacts the up-the-ramp samples by making them \"[jump](https://www.youtube.com/watch?v=SwYN7mTi6HM)\" from one group to another. We already noted this happening above\n", "[when we discussed saturation](#saturation) --- a pixel was suddenly pushed above the saturation limit and the `saturation` step flagged the pixel. However, some other jumps are not as dramatic, and the data after the jump might actually be as usable as data before the jump.\n", "\n", - "\n", "#### 4.7.1 Understanding jumps and the `jump` step\n", "\n", - "To exemplify the behavior of the jumps in up-the-ramp samples, let's look at an example. Consider the behavior of pixel index `(12,1000)` in integration `67`:" + "To exemplify the behavior of the jumps in up-the-ramp samples, let's look at an example. Consider the behavior of pixel index `(12,1000)` in integration `67`:\n" ] }, { @@ -2102,7 +2250,9 @@ "metadata": {}, "outputs": [], "source": [ - "jump_results = calwebb_detector1.jump_step.JumpStep.call(refpix_results, maximum_cores='all')" + "jump_results = calwebb_detector1.jump_step.JumpStep.call(\n", + " refpix_results, maximum_cores=\"all\"\n", + ")" ] }, { @@ -2132,22 +2282,39 @@ "column_index = 213\n", "row_index = 182\n", "\n", - "plt.title(f'Pixel index ({row_index}, {column_index})')\n", + "plt.title(f\"Pixel index ({row_index}, {column_index})\")\n", "\n", "group = np.arange(uncal_nis[0].data.shape[1])\n", "\n", - "plt.plot(group+1, uncal_nis[0].data[67, :, row_index, column_index], 'o-', \n", - " color='black', mfc='white', label='Integration 67'\n", - " )\n", - "plt.plot(group+1, uncal_nis[0].data[66, :, row_index, column_index], 'o-', \n", - " color='tomato', mfc='white', label='Integration 66', alpha=0.5\n", - " )\n", - "plt.plot(group+1, uncal_nis[0].data[68, :, row_index, column_index], 'o-', \n", - " color='cornflowerblue', mfc='white', label='Integration 68', alpha=0.5\n", - " )\n", - "\n", - "plt.xlabel('Group number', fontsize=16)\n", - "plt.ylabel('Counts', fontsize=16)\n", + "plt.plot(\n", + " group + 1,\n", + " uncal_nis[0].data[67, :, row_index, column_index],\n", + " \"o-\",\n", + " color=\"black\",\n", + " mfc=\"white\",\n", + " label=\"Integration 67\",\n", + ")\n", + "plt.plot(\n", + " group + 1,\n", + " uncal_nis[0].data[66, :, row_index, column_index],\n", + " \"o-\",\n", + " color=\"tomato\",\n", + " mfc=\"white\",\n", + " label=\"Integration 66\",\n", + " alpha=0.5,\n", + ")\n", + "plt.plot(\n", + " group + 1,\n", + " uncal_nis[0].data[68, :, row_index, column_index],\n", + " \"o-\",\n", + " color=\"cornflowerblue\",\n", + " mfc=\"white\",\n", + " label=\"Integration 68\",\n", + " alpha=0.5,\n", + ")\n", + "\n", + "plt.xlabel(\"Group number\", fontsize=16)\n", + "plt.ylabel(\"Counts\", fontsize=16)\n", "plt.legend()\n", "plt.show()" ] @@ -2157,7 +2324,7 @@ "id": "b0bde0f8-c2ad-4e3c-ac64-7805f8083ac6", "metadata": {}, "source": [ - "While the intercept of the different up-the-ramp samples is slightly different, the _slope_ (i.e., the count-rate) of it is fairly similar for integrations 66, 67 and 68. However, integration 67 shows a clear jump at group 4, likely from a cosmic ray. Let's take a look at what happened in this integration and group in the 2D spectrum:" + "While the intercept of the different up-the-ramp samples is slightly different, the _slope_ (i.e., the count-rate) of it is fairly similar for integrations 66, 67 and 68. However, integration 67 shows a clear jump at group 4, likely from a cosmic ray. Let's take a look at what happened in this integration and group in the 2D spectrum:\n" ] }, { @@ -2175,23 +2342,23 @@ "plt.subplot(1, 3, 1)\n", "im = plt.imshow(uncal_nis[0].data[i_integration, i_group, :, :])\n", "im.set_clim(-100, 1000)\n", - "plt.xlim(column_index-5, column_index+5)\n", - "plt.ylim(row_index-5, row_index+5)\n", - "plt.title('Integration 66, group 15')\n", + "plt.xlim(column_index - 5, column_index + 5)\n", + "plt.ylim(row_index - 5, row_index + 5)\n", + "plt.title(\"Integration 66, group 15\")\n", "\n", "plt.subplot(1, 3, 2)\n", - "im = plt.imshow(uncal_nis[0].data[i_integration+1, i_group, :, :])\n", + "im = plt.imshow(uncal_nis[0].data[i_integration + 1, i_group, :, :])\n", "im.set_clim(-100, 1000)\n", - "plt.xlim(column_index-5, column_index+5)\n", - "plt.ylim(row_index-5, row_index+5)\n", - "plt.title('Integration 67, group 15')\n", + "plt.xlim(column_index - 5, column_index + 5)\n", + "plt.ylim(row_index - 5, row_index + 5)\n", + "plt.title(\"Integration 67, group 15\")\n", "\n", "plt.subplot(1, 3, 3)\n", - "im = plt.imshow(uncal_nis[0].data[i_integration+2, i_group, :, :])\n", + "im = plt.imshow(uncal_nis[0].data[i_integration + 2, i_group, :, :])\n", "im.set_clim(-100, 1000)\n", - "plt.xlim(column_index-5, column_index+5)\n", - "plt.ylim(row_index-5, row_index+5)\n", - "plt.title('Integration 68, group 15')\n", + "plt.xlim(column_index - 5, column_index + 5)\n", + "plt.ylim(row_index - 5, row_index + 5)\n", + "plt.title(\"Integration 68, group 15\")\n", "plt.show()" ] }, @@ -2200,7 +2367,7 @@ "id": "d824c640-e5a7-4292-a5c3-d53bb1ec0ff1", "metadata": {}, "source": [ - "Ah! Clearly some cosmic ray hitting around pixel `(182, 213)`, with an area of about a pixel --- including pixel `(182, 213)`. Note that the `groupdq` doesn't show anything unusual so far:" + "Ah! Clearly some cosmic ray hitting around pixel `(182, 213)`, with an area of about a pixel --- including pixel `(182, 213)`. Note that the `groupdq` doesn't show anything unusual so far:\n" ] }, { @@ -2218,9 +2385,9 @@ "id": "206b1c74", "metadata": {}, "source": [ - "The JWST Calibration pipeline has an algorithm that aims to detect those jumps --- and is appropriately named the `jump` step. An important consideration when running the `jump` step is that one can use multiprocessing to run the step. This can offer dramatic speed improvements when running the step, in particular on large subarrays of data. The number of cores to use can be defined by the `maximum_cores` parameter, which can be an integer number or `all`, which will use all available cores. \n", + "The JWST Calibration pipeline has an algorithm that aims to detect those jumps --- and is appropriately named the `jump` step. An important consideration when running the `jump` step is that one can use multiprocessing to run the step. This can offer dramatic speed improvements when running the step, in particular on large subarrays of data. The number of cores to use can be defined by the `maximum_cores` parameter, which can be an integer number or `all`, which will use all available cores.\n", "\n", - "Let's run the step using all cores (this step does take some time ~4 mins):" + "Let's run the step using all cores (this step does take some time ~4 mins):\n" ] }, { @@ -2231,7 +2398,9 @@ "outputs": [], "source": [ "for i in range(nsegments):\n", - " uncal_nis[i] = calwebb_detector1.jump_step.JumpStep.call(uncal_nis[i], maximum_cores='all')" + " uncal_nis[i] = calwebb_detector1.jump_step.JumpStep.call(\n", + " uncal_nis[i], maximum_cores=\"all\"\n", + " )" ] }, { @@ -2239,7 +2408,7 @@ "id": "a0d989a3-62f1-4ab2-ba3d-73146bea5396", "metadata": {}, "source": [ - "It's not too obvious from the messages in the pipeline what happened, but the algorithm was used to _detect_ jumps, and these are added as new data-quality flags in the `groupdq`. Let's see what happened with the pixel identified by eye above:" + "It's not too obvious from the messages in the pipeline what happened, but the algorithm was used to _detect_ jumps, and these are added as new data-quality flags in the `groupdq`. Let's see what happened with the pixel identified by eye above:\n" ] }, { @@ -2267,7 +2436,7 @@ "id": "5b64bf86", "metadata": {}, "source": [ - "Aha! It changed. What does this mean? Let's repeat the trick we learned with the `saturation` step:" + "Aha! It changed. What does this mean? Let's repeat the trick we learned with the `saturation` step:\n" ] }, { @@ -2278,9 +2447,9 @@ "outputs": [], "source": [ "datamodels.dqflags.dqflags_to_mnemonics(\n", - " uncal_nis[0].groupdq[67, -1, row_index, column_index], \n", - " mnemonic_map=datamodels.dqflags.pixel\n", - " )" + " uncal_nis[0].groupdq[67, -1, row_index, column_index],\n", + " mnemonic_map=datamodels.dqflags.pixel,\n", + ")" ] }, { @@ -2288,11 +2457,11 @@ "id": "6a487a84", "metadata": {}, "source": [ - "Nice! We now have a flag that identifies when a jump detection happened. \n", + "Nice! We now have a flag that identifies when a jump detection happened.\n", "\n", "#### 4.7.2 Jump rates per integration\n", "\n", - "For fun, let's use the `groupdq` changes to figure out how many jumps happened per integration on this first segment of data by simple differencing with the products from the previous step, the `dark_current` step:" + "For fun, let's use the `groupdq` changes to figure out how many jumps happened per integration on this first segment of data by simple differencing with the products from the previous step, the `dark_current` step:\n" ] }, { @@ -2308,8 +2477,11 @@ "# Iterate through integrations counting how many pixels changed in all groups:\n", "for integration in range(uncal_nis[0].groupdq.shape[0]):\n", "\n", - " groupdq_difference = uncal_nis[0].groupdq[integration, :, :, :] - darkcurrent_results.groupdq[integration, :, :, :]\n", - " wherejumps = np.where(groupdq_difference != 0.)\n", + " groupdq_difference = (\n", + " uncal_nis[0].groupdq[integration, :, :, :]\n", + " - darkcurrent_results.groupdq[integration, :, :, :]\n", + " )\n", + " wherejumps = np.where(groupdq_difference != 0.0)\n", " njumps[integration] = len(wherejumps[0])" ] }, @@ -2318,7 +2490,7 @@ "id": "3b238f93-18ae-4453-9e8b-9962e87cb075", "metadata": {}, "source": [ - "Let's plot this:" + "Let's plot this:\n" ] }, { @@ -2331,10 +2503,10 @@ "integrations = np.arange(uncal_nis[0].groupdq.shape[0]) + 1\n", "\n", "plt.figure(figsize=(6, 4))\n", - "plt.title('Number of jumps on the first segment of data for NIS')\n", - "plt.plot(integrations, njumps, 'o-', color='black', mfc='white')\n", - "plt.xlabel('Integration', fontsize=16)\n", - "plt.ylabel('Number of jumps', fontsize=16)\n", + "plt.title(\"Number of jumps on the first segment of data for NIS\")\n", + "plt.plot(integrations, njumps, \"o-\", color=\"black\", mfc=\"white\")\n", + "plt.xlabel(\"Integration\", fontsize=16)\n", + "plt.ylabel(\"Number of jumps\", fontsize=16)\n", "plt.xlim(0.5, uncal_nis[0].groupdq.shape[0] + 0.5)\n", "plt.show()" ] @@ -2344,7 +2516,7 @@ "id": "e8b104d5-0b46-4ddb-a0ec-afa10717004f", "metadata": {}, "source": [ - "Very interesting! Per integration, it seems on the order of ~3,500 average jumps are detected. Each integration has (ngroups) x (number of pixels) = 70 x 32 x 2048 = 4587520 opportunities for jumps to appear, so this means an average rate of (detected events) / (total opportunities) = 0.07% per integration for this particular segment, detector and dataset." + "Very interesting! Per integration, it seems on the order of ~3,500 average jumps are detected. Each integration has (ngroups) x (number of pixels) = 70 x 32 x 2048 = 4587520 opportunities for jumps to appear, so this means an average rate of (detected events) / (total opportunities) = 0.07% per integration for this particular segment, detector and dataset.\n" ] }, { @@ -2352,7 +2524,7 @@ "id": "f50a6cf0-3433-40dc-90b7-b579ba115f10", "metadata": {}, "source": [ - "
Note on the effectiveness of the jump detection step: The jump detection step uses, by default, a two-point difference method that relies on appropriate knowledge of the read-noise of the detector. In some cases, this might be significantly off (or detector1 corrections might not be optimal as to leave significant detector effects) such that the algorithm might be shown to be too aggressive. Similarly, the algorithm relies on a decent amount of groups in the integration to work properly (larger than about 5). It is, thus, important to try different parameters to identify jumps in a given dataset and study their impact on the final products. One of the most important parameters is the rejection_threshold. The default value is 4, but TSO studies in the literature have sometimes opted for more conservative values (typically larger than 10). For this particular dataset, which has a large number of groups (70), the default value works well, but it might not be optimal nor be the best for other datasets." + "
Note on the effectiveness of the jump detection step: The jump detection step uses, by default, a two-point difference method that relies on appropriate knowledge of the read-noise of the detector. In some cases, this might be significantly off (or detector1 corrections might not be optimal as to leave significant detector effects) such that the algorithm might be shown to be too aggressive. Similarly, the algorithm relies on a decent amount of groups in the integration to work properly (larger than about 5). It is, thus, important to try different parameters to identify jumps in a given dataset and study their impact on the final products. One of the most important parameters is the rejection_threshold. The default value is 4, but TSO studies in the literature have sometimes opted for more conservative values (typically larger than 10). For this particular dataset, which has a large number of groups (70), the default value works well, but it might not be optimal nor be the best for other datasets.\n" ] }, { @@ -2360,7 +2532,7 @@ "id": "c3afc925-f0b2-49dc-80af-ec9ef301f5a8", "metadata": {}, "source": [ - "Before moving to the next step, we showcase one additional function from the `datamodels` which allows to save products to files --- the `save` function. This step is optional, if you want to use the `jump` step products for later use, uncomment the lines below:" + "Before moving to the next step, we showcase one additional function from the `datamodels` which allows to save products to files --- the `save` function. This step is optional, if you want to use the `jump` step products for later use, uncomment the lines below:\n" ] }, { @@ -2370,7 +2542,7 @@ "metadata": {}, "outputs": [], "source": [ - "# uncomment this line to run \n", + "# uncomment this line to run\n", "# if not os.path.exists('nis_jumpstep_seg001.fits'):\n", "# nsegments = 3\n", "# for i in range(nsegments):\n", @@ -2385,13 +2557,13 @@ "source": [ "### 4.8 Fitting ramps with the `ramp_fit` 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. The reason why this is not straightforward to do is because samples up-the-ramp are correlated. That is, because signal is accumulated up-the-ramp, group number 2 has a non-zero covariance with group number 1, and so on. \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. The reason why this is not straightforward to do is because samples up-the-ramp are correlated. That is, because signal is accumulated up-the-ramp, group number 2 has a non-zero covariance with group number 1, and so on.\n", "\n", "In addition, we will save the results to this step to a desired `output_dir` location and used in the next notebook where we'll generate some Lightcurvs.\n", "\n", "#### 4.8.1 Applying the `ramp_fit` step to JWST data\n", "\n", - "The JWST Calibration pipeline algorithm performs a sensible weighting of each group in order to account for that correlation when fitting a slope on the samples. Let's run this step, and save the products in files as we go, so we can use them for the next notebook. Note that as in the `jump` step, we can also run this step via multi-processing --- and we do just that below (if not ran already):" + "The JWST Calibration pipeline algorithm performs a sensible weighting of each group in order to account for that correlation when fitting a slope on the samples. Let's run this step, and save the products in files as we go, so we can use them for the next notebook. Note that as in the `jump` step, we can also run this step via multi-processing --- and we do just that below (if not ran already):\n" ] }, { @@ -2409,11 +2581,8 @@ "\n", "for i in range(nsegments):\n", " uncal_nis[i] = calwebb_detector1.ramp_fit_step.RampFitStep.call(\n", - " uncal_nis[i], \n", - " maximum_cores='all', \n", - " save_results=True, \n", - " output_dir=output_dir\n", - " )" + " uncal_nis[i], maximum_cores=\"all\", save_results=True, output_dir=output_dir\n", + " )" ] }, { @@ -2421,7 +2590,7 @@ "id": "96724249", "metadata": {}, "source": [ - "All right, note the products of this step for TSO's are actually a list:" + "All right, note the products of this step for TSO's are actually a list:\n" ] }, { @@ -2439,7 +2608,7 @@ "id": "159fe870", "metadata": {}, "source": [ - "The data associated with the zeroth element of this list (`ramps_nis1[0][0].data`) has dimensions equal to the size of the frames (rows and columns). The first element (`ramps_nis1[0][1].data`), has three dimensions, the same as the zeroth but for each integration. We usually refer to this latter product as the `rateints` product --- i.e., the rates per integration:" + "The data associated with the zeroth element of this list (`ramps_nis1[0][0].data`) has dimensions equal to the size of the frames (rows and columns). The first element (`ramps_nis1[0][1].data`), has three dimensions, the same as the zeroth but for each integration. We usually refer to this latter product as the `rateints` product --- i.e., the rates per integration:\n" ] }, { @@ -2467,7 +2636,7 @@ "id": "34ae6b9d", "metadata": {}, "source": [ - "To familiarize ourselves with these products, let's plot the rates of the 10th integration for NIS:" + "To familiarize ourselves with these products, let's plot the rates of the 10th integration for NIS:\n" ] }, { @@ -2478,10 +2647,10 @@ "outputs": [], "source": [ "plt.figure(figsize=(12, 3))\n", - "plt.title('NIS data; rates for integration 10')\n", + "plt.title(\"NIS data; rates for integration 10\")\n", "im = plt.imshow(uncal_nis[0][1].data[10, :, :])\n", "im.set_clim(-1, 10)\n", - "plt.colorbar(label='Counts/s')\n", + "plt.colorbar(label=\"Counts/s\")\n", "plt.show()" ] }, @@ -2490,7 +2659,7 @@ "id": "3c7fcc50-a114-4c2d-a487-a9f5e493b912", "metadata": {}, "source": [ - "In case you were unsure of the units in the colorbar, you can double-check them through the `datamodels` themselves:" + "In case you were unsure of the units in the colorbar, you can double-check them through the `datamodels` themselves:\n" ] }, { @@ -2500,7 +2669,7 @@ "metadata": {}, "outputs": [], "source": [ - "uncal_nis[0][1].search('unit')" + "uncal_nis[0][1].search(\"unit\")" ] }, { @@ -2508,7 +2677,7 @@ "id": "8998289d", "metadata": {}, "source": [ - "These rates look very pretty, lets check the first element results for the 10th inegration." + "These rates look very pretty, lets check the first element results for the 10th inegration.\n" ] }, { @@ -2519,10 +2688,10 @@ "outputs": [], "source": [ "plt.figure(figsize=(12, 3))\n", - "plt.title('NIS data; rates for integration 10')\n", + "plt.title(\"NIS data; rates for integration 10\")\n", "im = plt.imshow(uncal_nis[0][1].data[10, :, :])\n", "im.set_clim(-1, 30)\n", - "plt.colorbar(label='Counts/s')\n", + "plt.colorbar(label=\"Counts/s\")\n", "plt.show()" ] }, @@ -2531,7 +2700,7 @@ "id": "dc22644d", "metadata": {}, "source": [ - "These rates look _very_ good as well. " + "These rates look _very_ good as well.\n" ] }, { @@ -2541,7 +2710,7 @@ "source": [ "## 5. Final words \n", "\n", - "This completes this notebook where we have reduced and calibrated NIRISS/SOSS data of WASP-39b from program 1366 using STAGE1 of the JWST pipeline. In the next notebook, `02_niriss_soss_spec2_generate_lightcurves`, we will use the calibrated data products to extract the spectra of WASP-39b and generate some lightcurves performing similiar steps to what is done in STAGE 2 of the JWST pipeline. I would like to thank the JWST NIRISS team, especially Néstor Espinoza and Aarynn Carter for their feedback and support For this particular effort of writing these NIRISS/SOSS demonstration notebooks." + "This completes this notebook where we have reduced and calibrated NIRISS/SOSS data of WASP-39b from program 1366 using STAGE1 of the JWST pipeline. In the next notebook, `02_niriss_soss_spec2_generate_lightcurves`, we will use the calibrated data products to extract the spectra of WASP-39b and generate some lightcurves performing similiar steps to what is done in STAGE 2 of the JWST pipeline. I would like to thank the JWST NIRISS team, especially Néstor Espinoza and Aarynn Carter for their feedback and support For this particular effort of writing these NIRISS/SOSS demonstration notebooks.\n" ] } ],