{ "cells": [ { "cell_type": "markdown", "id": "3a0cc965", "metadata": {}, "source": [ "\n", "# NIRSpec MOS pipeline processing" ] }, { "cell_type": "markdown", "id": "df4f60f7", "metadata": {}, "source": [ "**Author**: James Muzerolle\n", "\n", "Plotting function originally developed by Bryan Hilbert\n", "\n", "**Latest Update**: 28 October 2021" ] }, { "cell_type": "markdown", "id": "edbd4add", "metadata": {}, "source": [ "## Table of Contents\n", "* [Introduction](#intro)\n", " * [Overview](#overview)\n", " * [Simulated data](#sims)\n", "* [Imports](#imports)\n", "* [Convenience functions](#func)\n", "* [Pipeline processing flag](#flag)\n", "* [Input simulations](#inputs)\n", "* [Association files and metadata](#associations)\n", "* [Run the calwebb_spec2 pipeline](#runspec2)\n", "* [Run the calwebb_spec3 pipeline](#runspec3)\n", "* [Master background subtraction](#mbs)" ] }, { "cell_type": "markdown", "id": "a7eb5dc4", "metadata": {}, "source": [ "## Introduction " ] }, { "cell_type": "markdown", "id": "58e6e77d", "metadata": {}, "source": [ "In this notebook, we will explore the stage 2 and 3 pipelines for NIRSpec MOS data. Here, we will focus on the mechanics of processing \"real\" example data, including how to use associations for exposure specification and multi-exposure combination, the role of metadata, and what the primary data products at each stage look like. We will also see examples of how to interact and work with data models and metadata.\n", "\n", "We are using pipeline version 1.2.3 for all data processing in this notebook. Most of the processing runs shown here use the default reference files from the Calibration Reference Data System (CRDS). Please note that pipeline software development is a continuous process, so results in some cases may be slightly different if using a different version. There are also a few known issues with some of the pipeline steps in this build that are expected to be fixed in the near future, though these do not significantly effect the products you will see here. Finally, some of the calibration reference files used by individual pipeline steps in the current CRDS context are placeholders, as they require calibrations that can only be taken in flight; this means that the absolute flux values and formal uncertainties in the products shown here are not physical and should not be taken literally." ] }, { "cell_type": "markdown", "id": "54c717ff", "metadata": {}, "source": [ "### Simulated data " ] }, { "cell_type": "markdown", "id": "f775a438", "metadata": {}, "source": [ "During the webinar, we used simulated NIRSpec exposures generated by the ESA Instrument Performance Simulator (IPS). The simulations were created using an input scene consisting of 60 point sources within the MSA field of view, each described by an emission-line galaxy template spectrum at different redshifts. The observation included a series of three nodded exposures, in which each source is moved to a different shutter within its three-shutter slitlet. The instrument was configured with two different disperser settings: G235M+F170LP and G395M+F290LP. The nodded exposures are numbered for convenience in these example file names according to the ordering of the positions in the pattern. These file names are not indicative of actual data product file names you will see in the archive.\n", "\n", "There are a number of caveats to be aware of regarding these simulated data. 1) The IPS does not include a full treatment of all of the effects corrected by the stage 2 pipeline, particularly some of the throughput components. As with the above caveat regarding reference files, the simulations shown here should not be used for any analyses of flux information. 2) The simulated PSF is truncated in order to save on computation time, which results in an artifical drop in apparent count rate in the PSF wings in some cases. 3) Some of the error arrays in the rate files are not populated, as those are normally generated during stage 1 of the pipeline. In addition to the reference file issue noted above, this means the final combined error values should be ignored. 4) Spacecraft pointing-related information has to be added by-hand to the headers before the simulated data can be processed by the pipeline. These keywords are used by the stage 3 pipeline when combining exposures, in order to align the spectral traces. In addition, the MSA metafile, used by the pipeline to determine open shutter and source information, also had to be created by-hand. Because this has to be a manual process and may be subject to small errors, the quality of the combined products here should not be taken as indicative of actual in-flight performance.\n", "\n", "***Because these simulations are not public yet, we cannot share them in the main repository. We will update the notebook with an appropriate dataset in the near future. For now, if you have a simulation, you will have to input it in the appropriate cell.***" ] }, { "cell_type": "markdown", "id": "1ab1a43a", "metadata": {}, "source": [ "## Imports " ] }, { "cell_type": "markdown", "id": "bcf18c69", "metadata": {}, "source": [ "Import packages necessary for this notebook" ] }, { "cell_type": "code", "execution_count": null, "id": "c9d45c63", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "import glob\n", "\n", "import os\n", "import zipfile\n", "import urllib.request\n", "\n", "import json\n", "\n", "from shutil import copyfile\n", "\n", "from astropy.io import fits\n", "from astropy.utils.data import download_file\n", "import astropy.units as u\n", "from astropy import wcs\n", "from astropy.wcs import WCS\n", "from astropy.visualization import ImageNormalize, ManualInterval, LogStretch, LinearStretch, AsinhStretch" ] }, { "cell_type": "markdown", "id": "39940c7b", "metadata": {}, "source": [ "Set up matplotlib for plotting" ] }, { "cell_type": "code", "execution_count": null, "id": "51f52309", "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import matplotlib as mpl\n", "\n", "# Use this version for non-interactive plots (easier scrolling of the notebook)\n", "%matplotlib inline\n", "\n", "# Use this version (outside of Jupyter Lab) if you want interactive plots\n", "#%matplotlib notebook\n", "\n", "# These gymnastics are needed to make the sizes of the figures\n", "# be the same in both the inline and notebook versions\n", "%config InlineBackend.print_figure_kwargs = {'bbox_inches': None}\n", "\n", "mpl.rcParams['savefig.dpi'] = 80\n", "mpl.rcParams['figure.dpi'] = 80\n", "\n", "plt.rcParams.update({'font.size': 18})" ] }, { "cell_type": "markdown", "id": "af2bb9f0", "metadata": {}, "source": [ "Import JWST pipeline modules" ] }, { "cell_type": "code", "execution_count": null, "id": "63872971", "metadata": {}, "outputs": [], "source": [ "# The calwebb_spec and spec3 pipelines\n", "from jwst.pipeline import Spec2Pipeline\n", "from jwst.pipeline import Spec3Pipeline\n", "\n", "# individual steps\n", "from jwst.assign_wcs import AssignWcsStep\n", "from jwst.assign_wcs import nirspec\n", "from jwst.background import BackgroundStep\n", "from jwst.imprint import ImprintStep\n", "from jwst.msaflagopen import MSAFlagOpenStep\n", "from jwst.extract_2d import Extract2dStep\n", "from jwst.srctype import SourceTypeStep\n", "from jwst.wavecorr import WavecorrStep\n", "from jwst.flatfield import FlatFieldStep\n", "from jwst.pathloss import PathLossStep\n", "from jwst.photom import PhotomStep\n", "from jwst.cube_build import CubeBuildStep\n", "from jwst.extract_1d import Extract1dStep\n", "from jwst.combine_1d import Combine1dStep\n", "\n", "# data models\n", "from jwst import datamodels\n", "\n", "# associations\n", "from jwst.associations import asn_from_list\n", "from jwst.associations.lib.rules_level3_base import DMS_Level3_Base" ] }, { "cell_type": "markdown", "id": "170a090f", "metadata": {}, "source": [ "## Define convenience functions and parameters \n", "\n", "**Update the link and file names to point at your simulation. This cell will fail if you do not input a valid path.**" ] }, { "cell_type": "code", "execution_count": null, "id": "470200ee", "metadata": {}, "outputs": [], "source": [ "# All files created by the steps in this notebook have been pre-computed and cached on Box\n", "\n", "# first specify a desired local directory in which to place the downloaded data, as well as any offline processing you choose to run\n", "output_dir = 'nirspec_files/'\n", "if not os.path.exists(output_dir):\n", " os.makedirs(output_dir)\n", "\n", "# Update the link to point at your simulation\n", "ziplink = 'MySimulationURL'\n", "zipfilename = 'mydata.zip'\n", "if not os.path.isfile(os.path.join(output_dir, zipfilename)):\n", " print('Downloading {}...'.format(zipfilename))\n", " demo_file = download_file(ziplink, cache=True)\n", " # Make a symbolic link using a local name for convenience\n", " os.symlink(demo_file, os.path.join(output_dir, zipfilename))\n", "else:\n", " print('{} already exists, skipping download...'.format(zipfilename))\n", "\n", "# unzip\n", "zf = zipfile.ZipFile(output_dir+zipfilename, 'r')\n", "print(output_dir)\n", "zf.extractall(output_dir)" ] }, { "cell_type": "code", "execution_count": null, "id": "fbc8efc2", "metadata": {}, "outputs": [], "source": [ "def show_image(data_2d, vmin, vmax, xsize=19, ysize=10, title=None, aspect=1, scale='log', units='MJy/sr'):\n", " \"\"\"Function to generate a 2D, log-scaled image of the data\n", " \n", " Parameters\n", " ----------\n", " data_2d : numpy.ndarray\n", " 2D image to be displayed\n", " \n", " vmin : float\n", " Minimum signal value to use for scaling\n", " \n", " vmax : float\n", " Maximum signal value to use for scaling\n", " \n", " title : str\n", " String to use for the plot title\n", " \n", " scale : str\n", " Specify scaling of the image. Can be 'log' or 'linear'\n", " \n", " units : str\n", " Units of the data. Used for the annotation in the\n", " color bar\n", " \"\"\"\n", " if scale == 'log':\n", " norm = ImageNormalize(data_2d, interval=ManualInterval(vmin=vmin, vmax=vmax),\n", " stretch=LogStretch())\n", " elif scale == 'linear':\n", " norm = ImageNormalize(data_2d, interval=ManualInterval(vmin=vmin, vmax=vmax),\n", " stretch=LinearStretch())\n", " elif scale == 'Asinh':\n", " norm = ImageNormalize(data_2d, interval=ManualInterval(vmin=vmin, vmax=vmax),\n", " stretch=AsinhStretch())\n", " fig = plt.figure(figsize=(xsize, ysize))\n", " ax = fig.add_subplot(1, 1, 1)\n", " im = ax.imshow(data_2d, origin='lower', norm=norm, aspect=aspect, cmap='gist_earth')\n", "\n", " if (units != 'none'):\n", " fig.colorbar(im, label=units)\n", " plt.xlabel('Pixel column')\n", " plt.ylabel('Pixel row')\n", " if title:\n", " plt.title(title)" ] }, { "cell_type": "markdown", "id": "ac500fc0", "metadata": {}, "source": [ "## Pipeline processing flag " ] }, { "cell_type": "markdown", "id": "02552e8c", "metadata": {}, "source": [ "The pipeline and individual steps take too long to run in real time for this demo, so all products shown here have been pre-computed, and the actual pipeline calls will be skipped. Change the following flag to True if you want to run everything offline." ] }, { "cell_type": "code", "execution_count": null, "id": "c2f237b9", "metadata": {}, "outputs": [], "source": [ "runflag = False" ] }, { "cell_type": "markdown", "id": "d08ce68d", "metadata": {}, "source": [ "## Input simulations " ] }, { "cell_type": "markdown", "id": "185481c6", "metadata": {}, "source": [ "At the webinar, we used simulated NIRSpec exposures of point sources observed with the MSA. Because the simulator generates count rate maps, equivalent to level 2a data products, we have to skip stage 1 of the pipeline and instead start the processing with the calwebb_spec2 pipeline. First, let's take a look at a few of the level 2a images to get familiarized with the inputs. The observation consists of 3 nodded exposures, with each source moving to a different shutter in its 3-shutter slitlet.\n", "\n", "**Update the file names to your simulated dataset here and in following cells.**" ] }, { "cell_type": "code", "execution_count": null, "id": "b0a03990", "metadata": {}, "outputs": [], "source": [ "# get the data model of dither position 1:\n", "ratefile1 = output_dir+'mydata_nrs1_rate.fits' # for the NRS1 detector\n", "dither1 = datamodels.open(ratefile1)\n", "ratefile2 = output_dir+'mydata_nrs2_rate.fits' # for the NRS2 detector\n", "dither2 = datamodels.open(ratefile2)\n", "\n", "# get the pixel data (the SCI extension of the fits file)\n", "ratesci1 = dither1.data\n", "ratesci2 = dither2.data\n", "\n", "# display the images\n", "show_image(ratesci1, 0, 1., xsize=19, ysize=19, units='DN/s', title='NRS1')\n", "show_image(ratesci2, 0, 1., xsize=19, ysize=19, units='DN/s', title='NRS2')\n", "\n", "# zoom in to show details of one slitlet\n", "show_image(ratesci1[1550:1650, 850:1050], 0, 0.1, xsize=19, ysize=15, units='none', scale='Asinh', title='zoom of NRS1')" ] }, { "cell_type": "markdown", "id": "b701e4cf", "metadata": {}, "source": [ "## Association files and metadata " ] }, { "cell_type": "markdown", "id": "371edc36", "metadata": {}, "source": [ "One purpose of the nodded exposures in this case is to enable background subtraction. The spec2 pipeline is set up to handle this by using association files that list the exposures to be used as backgrounds for each input exposure. The background elements can also be specified manually as inputs into the \"background\" step, but we'll be using association files in this example." ] }, { "cell_type": "code", "execution_count": null, "id": "01b45229", "metadata": {}, "outputs": [], "source": [ "# show the contents of one of the association files\n", "asn_file = output_dir+\"l2_g235m_asn.json\"\n", "with open(asn_file) as f_obj:\n", " asn_data = json.load(f_obj)\n", "asn_data" ] }, { "cell_type": "markdown", "id": "d237f24c", "metadata": {}, "source": [ "One unique aspect of MOS data processing is that it requires additional metadata, such as shutter locations and source positions. The pipeline gets this information from an MSA \"metafile\" that is automatically generated, along with the raw \"level 1b\" files, using metadata taken from the PPSDB. The metadata is originally populated when the observation is created with the MSA planning tool in APT. The metafile is a fits file containing two binary fits tables, one with MSA shutter information for each slitlet and nod, and the other with various details about each source observed. Let's take a look at the metafile for this simulation:" ] }, { "cell_type": "code", "execution_count": null, "id": "90404df8", "metadata": {}, "outputs": [], "source": [ "# open the MSA metafile\n", "metafile = output_dir+'pointing_summary_n_metafile_msa.fits'\n", "hdul = fits.open(metafile)\n", "hdul.info()" ] }, { "cell_type": "code", "execution_count": null, "id": "4c958f5e", "metadata": {}, "outputs": [], "source": [ "# columns of the \"SHUTTER INFO\" table\n", "hdul['SHUTTER_INFO'].columns" ] }, { "cell_type": "markdown", "id": "6432be8f", "metadata": {}, "source": [ "In the \"SHUTTER_INFO\" table, the number of rows per slitlet is N_s x N_n, where N_s is the number of open shutters in each slitlet and N_n is the number of nods. In this case, N_s = 3 and N_n = 3, so there are 9 rows per slitlet. These are needed by the pipeline to determine which slitlet shutter contains the source in each nodded exposure." ] }, { "cell_type": "code", "execution_count": null, "id": "a4910466", "metadata": {}, "outputs": [], "source": [ "# print table entries for the first two slitlets\n", "print(hdul['SHUTTER_INFO'].data[:18])" ] }, { "cell_type": "code", "execution_count": null, "id": "ecf31faf", "metadata": {}, "outputs": [], "source": [ "# columns of the \"SOURCE_INFO\" table\n", "hdul['SOURCE_INFO'].columns" ] }, { "cell_type": "markdown", "id": "b580df75", "metadata": {}, "source": [ "The SOURCE_INFO table contains one row per source that was observed with an MSA configuration. The most important information here are the source RA & Dec coordinates, and the stellarity parameter, which is used to determine whether the source should be processed as point or extended." ] }, { "cell_type": "code", "execution_count": null, "id": "034b3836", "metadata": {}, "outputs": [], "source": [ "# print table entries for the first 10 sources\n", "print(hdul['SOURCE_INFO'].data[:10])" ] }, { "cell_type": "markdown", "id": "9b0a6987", "metadata": {}, "source": [ "## Run the spec2 pipeline " ] }, { "cell_type": "code", "execution_count": null, "id": "d218ae96", "metadata": {}, "outputs": [], "source": [ "# run the calwebb_spec2 pipeline using an association file as input\n", "\n", "if runflag:\n", " spec2 = Spec2Pipeline()\n", " spec2.save_results = True\n", " spec2.output_dir = output_dir\n", " result = spec2(asn_file)" ] }, { "cell_type": "code", "execution_count": null, "id": "89f6b45e", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# take a look at the results - open the level 2b files\n", "\n", "callist = [f for f in glob.glob(output_dir+\"*g235m_nrs?_cal.fits\")]\n", "callist.sort()\n", "for calfile in callist:\n", " print(calfile)\n", " cal = datamodels.open(calfile) # this contains the calibrated unrectified 2D spectra\n", " root = calfile[:-9]\n", " s2d = datamodels.open(root+'_s2d.fits') # this contains the calibrated *rectified* 2D spectra\n", " x1d = datamodels.open(root+'_x1d.fits') # this contains the aperture-extracted 1D spectra\n", " \n", " for i, slit in enumerate(cal.slits):\n", "\n", " if (slit.name == '31'): # change this, or comment out, to see other slits\n", " print(slit.name)\n", " \n", " calsci = slit.data # contains the pixel data from the cal file (SCI extension)\n", " s2dsci = s2d.slits[i].data # contains the pixel data from the s2d file\n", " \n", " # determine the wavelength scale of the s2d data for plotting purposes\n", " # get the data model WCS object\n", " wcsobj = s2d.slits[i].meta.wcs\n", " y, x = np.mgrid[:s2dsci.shape[0], : s2dsci.shape[1]] # grid of pixel x,y indices\n", " det2sky = wcsobj.get_transform('detector', 'world') # the coordinate transform from detector space (pixels) to sky (RA, DEC in degrees)\n", " ra, dec, s2dwave = det2sky(x, y) # RA, Dec, wavelength (microns) for each pixel\n", " s2dwaves = s2dwave[0, :] # only need a single row of values since this is the rectified spectrum\n", " xtint = np.arange(100, s2dsci.shape[1], 100)\n", " xtlab = np.round(s2dwaves[xtint], 2) # wavelength labels for the x-axis\n", " \n", " # get wavelength & flux from the x1d data model\n", " x1dwave = x1d.spec[i].spec_table.WAVELENGTH\n", " x1dflux = x1d.spec[i].spec_table.FLUX\n", " \n", " # plot the unrectified calibrated 2D spectrum\n", " show_image(calsci, -0.01, 0.01, aspect=5., scale='linear', units='MJy')\n", " \n", " # plot the rectified 2D spectrum\n", " show_image(s2dsci, -0.01, 0.01, aspect=5., scale='linear', units='MJy')\n", " plt.xticks(xtint, xtlab)\n", " plt.xlabel('wavelength (microns)')\n", " \n", " # plot the 1D extracted spectrum\n", " fig = plt.figure(figsize=(19, 8))\n", " plt.plot(x1dwave, x1dflux)\n", " plt.xlabel('wavelength (microns)')\n", " plt.ylabel('flux (Jy)')\n", " plt.show()" ] }, { "cell_type": "markdown", "id": "0cda16e9", "metadata": {}, "source": [ "Now we'll process the G395M data, using the same methodology (to be used later)." ] }, { "cell_type": "code", "execution_count": null, "id": "1eaf7c00", "metadata": {}, "outputs": [], "source": [ "# association file\n", "asn_file = output_dir+\"l2_g395m_asn.json\"" ] }, { "cell_type": "code", "execution_count": null, "id": "8265cd99", "metadata": {}, "outputs": [], "source": [ "# run the calwebb_spec2 pipeline using association files as inputs\n", "\n", "if runflag:\n", " spec2 = Spec2Pipeline()\n", " spec2.save_results = True\n", " spec2.output_dir = output_dir\n", " result = spec2(asn_file)" ] }, { "cell_type": "markdown", "id": "48757537", "metadata": {}, "source": [ "## Run the calwebb_spec3 pipeline " ] }, { "cell_type": "markdown", "id": "d42ad9c8", "metadata": {}, "source": [ "For an observation that contains multiple exposures of the same sources, such as the nod set in this example, this pipeline will combine all of the exposures into a single output product. Both a rectified 2D and extracted 1D spectrum are generated. This process also includes an outlier detection step that compares the stack of values at each resampled pixel and flags outliers based on noise threshold parameters; the flagged outliers are not included in the spectral combination. However, we will skip that step in this example because of the aforementioned limitations with the noise propagation, as well as the fact that the simulations do not include outliers." ] }, { "cell_type": "code", "execution_count": null, "id": "cbf1911f", "metadata": {}, "outputs": [], "source": [ "# create an association file with all the level 2b cal.fits products generated above\n", "\n", "if runflag:\n", " outfilename = 'nirspec_mossim_g235m_{source_id}_combined.fits'\n", " calfilelist = [os.path.basename(f) for f in glob.glob(output_dir+'nirspec_mossim_d?_g235m_nrs?_cal.fits')]\n", " calfilelist.sort()\n", " \n", " asn = asn_from_list.asn_from_list(calfilelist, product_name=outfilename)\n", "\n", " asn_file_l3 = output_dir+'l3_g235m_asn.json'\n", " with open(asn_file_l3, 'w') as outfile:\n", " name, serialized = asn.dump(format='json')\n", " outfile.write(serialized)" ] }, { "cell_type": "code", "execution_count": null, "id": "b34b753c", "metadata": {}, "outputs": [], "source": [ "# run the calwebb_spec3 pipeline using the association file as input\n", "\n", "if runflag:\n", " spec3 = Spec3Pipeline()\n", " spec3.save_results = True\n", " spec3.output_dir = output_dir\n", " spec3.outlier_detection.skip = True\n", " result = spec3(asn_file_l3)" ] }, { "cell_type": "code", "execution_count": null, "id": "85edc90f", "metadata": {}, "outputs": [], "source": [ "# display the level 3 combined products\n", "\n", "s2d3list = [f for f in sorted(glob.glob(output_dir+'*g235m*_combined_s2d.fits'))]\n", "x1d3list = [f for f in sorted(glob.glob(output_dir+'*g235m*_combined_x1d.fits'))]\n", "\n", "for i, file in enumerate(s2d3list):\n", " s2d3 = datamodels.open(file)\n", " \n", " if (s2d3.name == '31'): # show just one of the files for brevity's sake; change or comment out to see others\n", " print(file)\n", " print(s2d3.name)\n", " s2d3sci = s2d3.data\n", " print(x1d3list[i])\n", " x1d3 = datamodels.open(x1d3list[i])\n", " x1dwave = x1d3.spec[0].spec_table.WAVELENGTH\n", " x1dflux = x1d3.spec[0].spec_table.FLUX\n", "\n", " # get the wavelength scale\n", " wcsobj = s2d3.meta.wcs\n", " y, x = np.mgrid[:s2d3sci.shape[0], : s2d3sci.shape[1]]\n", " det2sky = wcsobj.get_transform('detector', 'world')\n", " s2d3ra, s2d3dec, s2d3wave = det2sky(x, y)\n", " s2d3waves = s2d3wave[0, :]\n", " xtint = np.arange(100, s2d3sci.shape[1], 100)\n", " xtlab = np.round(s2d3waves[xtint], 2)\n", "\n", " # show the combined s2d image\n", " show_image(s2d3sci, -0.01, 0.01, aspect=5., scale='linear', units='MJy')\n", " plt.xticks(xtint, xtlab)\n", " plt.xlabel('wavelength (microns)')\n", "\n", " # plot extracted 1D spectrum\n", " fig = plt.figure(figsize=(19, 8))\n", " plt.plot(x1dwave, x1dflux)\n", " plt.xlabel('wavelength (microns)')\n", " plt.ylabel('flux (Jy)')\n", " plt.ylim(-5.e4, 6.e5)\n", " plt.show()" ] }, { "cell_type": "markdown", "id": "0f6580d2", "metadata": {}, "source": [ "Now let's combine the G395M data." ] }, { "cell_type": "code", "execution_count": null, "id": "f22b7519", "metadata": {}, "outputs": [], "source": [ "# create an association file with all the level 2b cal.fits products generated above\n", "\n", "if runflag:\n", " outfilename = 'nirspec_mossim_g395m_{source_id}_combined.fits'\n", " calfilelist = [os.path.basename(f) for f in glob.glob(output_dir+'nirspec_mossim_d?_g395m_nrs?_cal.fits')]\n", " calfilelist.sort()\n", " \n", " asn = asn_from_list.asn_from_list(calfilelist, product_name=outfilename)\n", "\n", " asn_file_l3 = output_dir+'l3_g395m_asn.json'\n", " with open(asn_file_l3, 'w') as outfile:\n", " name, serialized = asn.dump(format='json')\n", " outfile.write(serialized)" ] }, { "cell_type": "code", "execution_count": null, "id": "d7a1abc2", "metadata": {}, "outputs": [], "source": [ "# run the calwebb_spec3 pipeline using the association file as input\n", "\n", "if runflag:\n", " spec3 = Spec3Pipeline()\n", " spec3.save_results = True\n", " spec3.output_dir = output_dir\n", " spec3.outlier_detection.skip = True # skip this step for now, because the simulations do not include noise\n", " result = spec3(asn_file_l3)" ] }, { "cell_type": "code", "execution_count": null, "id": "9bf368f3", "metadata": {}, "outputs": [], "source": [ "# display the level 3 combined products\n", "\n", "s2d3list = [f for f in sorted(glob.glob(output_dir+'*g395m*_combined_s2d.fits'))]\n", "x1d3list = [f for f in sorted(glob.glob(output_dir+'*g395m*_combined_x1d.fits'))]\n", "\n", "for i, file in enumerate(s2d3list):\n", " s2d3 = datamodels.open(file)\n", " \n", " if (s2d3.name == '31'): # show just one of the files for brevity's sake; change or comment out to see others\n", " print(file)\n", " print(s2d3.name)\n", " s2d3sci = s2d3.data\n", " print(x1d3list[i])\n", " x1d3 = datamodels.open(x1d3list[i])\n", " x1dwave = x1d3.spec[0].spec_table.WAVELENGTH\n", " x1dflux = x1d3.spec[0].spec_table.FLUX\n", "\n", " # get the wavelength scale\n", " wcsobj = s2d3.meta.wcs\n", " y, x = np.mgrid[:s2d3sci.shape[0], : s2d3sci.shape[1]]\n", " det2sky = wcsobj.get_transform('detector', 'world')\n", " s2d3ra, s2d3dec, s2d3wave = det2sky(x, y)\n", " s2d3waves = s2d3wave[0, :]\n", " xtint = np.arange(100, s2d3sci.shape[1], 100)\n", " xtlab = np.round(s2d3waves[xtint], 2)\n", "\n", " # show the combined s2d image \n", " show_image(s2d3sci, -0.01, 0.01, aspect=5., scale='linear', units='MJy')\n", " plt.xticks(xtint, xtlab)\n", " plt.xlabel('wavelength (microns)')\n", "\n", " # plot extracted 1D spectrum\n", " fig = plt.figure(figsize=(19, 8))\n", " plt.plot(x1dwave, x1dflux)\n", " plt.xlabel('wavelength (microns)')\n", " plt.ylabel('flux (Jy)')\n", " plt.show()" ] }, { "cell_type": "markdown", "id": "e3cb6316", "metadata": {}, "source": [ "Now that we have level 3 spectra for both gratings, these can be combined at the 1D level." ] }, { "cell_type": "code", "execution_count": null, "id": "dee76ec6", "metadata": {}, "outputs": [], "source": [ "# association file with 1D spectra per source per grating\n", "asn_file = output_dir+'nirspec_mossim_all_s02140.json'\n", "\n", "if runflag:\n", " step = Combine1dStep()\n", " step.save_results = True\n", " step.output_dir = output_dir\n", " result = step(asn_file)" ] }, { "cell_type": "code", "execution_count": null, "id": "f410e194", "metadata": {}, "outputs": [], "source": [ "x1dcomb = datamodels.open(output_dir+'nirspec_mossim_all_s02140_combine1dstep.fits')\n", "x1dcombwave = x1dcomb.spec[0].spec_table.WAVELENGTH\n", "x1dcombflux = x1dcomb.spec[0].spec_table.FLUX\n", "\n", "# plot combined 1D spectrum\n", "fig = plt.figure(figsize=(19, 8))\n", "plt.plot(x1dcombwave, x1dcombflux)\n", "plt.xlabel('wavelength (microns)')\n", "plt.ylabel('flux (Jy)')\n", "plt.ylim(-5.e4, 2.e5)\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "e48f3fee", "metadata": {}, "source": [ "## Master background subtraction using MSA background shutters " ] }, { "cell_type": "markdown", "id": "254c977c", "metadata": {}, "source": [ "We'll now explore one example of the master background subraction methodology for MOS data. There are multiple ways of creating a master background; the most common will likely be programs that have designed specific background shutters for that purpose. The simulations we are using here do not include such dedicated backgrounds, but some of the targets have minimal S/N and hence are dominated by the background emission, so we will use those as a proxy after some minipulation of the MSA metafile." ] }, { "cell_type": "markdown", "id": "1876ea46", "metadata": {}, "source": [ "For comparison purposes, we will first reprocess one of the exposures without doing any background subtraction.\n", "To simplify for this demonstration, we will use just one of the sources observed in the exposure (source ID 2140, slitlet ID 31), isolating it by editing the MSA metafile appropriately to include only one row in the SHUTTER_INFO table." ] }, { "cell_type": "code", "execution_count": null, "id": "cf339397", "metadata": {}, "outputs": [], "source": [ "# edit the metafile to include only a single row corresponding to the source and exposure we want to use\n", "\n", "if runflag:\n", "\n", " metafile = output_dir+'pointing_summary_n_metafile_msa.fits'\n", " hdul = fits.open(metafile)\n", "\n", " # pick the row in the original metafile SHUTTER_INFO table with a reasonably bright source in the appropriate shutter for this particular nod position\n", " indices = [255]\n", "\n", " # find the matching values for each of the columns\n", " slitlets = hdul['SHUTTER_INFO'].data['SLITLET_ID'][indices]\n", " metaids = hdul['SHUTTER_INFO'].data['MSA_METADATA_ID'][indices]\n", " quads = hdul['SHUTTER_INFO'].data['SHUTTER_QUADRANT'][indices]\n", " rows = hdul['SHUTTER_INFO'].data['SHUTTER_ROW'][indices]\n", " cols = hdul['SHUTTER_INFO'].data['SHUTTER_COLUMN'][indices]\n", " sources = hdul['SHUTTER_INFO'].data['SOURCE_ID'][indices]\n", " bkgd = hdul['SHUTTER_INFO'].data['BACKGROUND'][indices]\n", " state = hdul['SHUTTER_INFO'].data['SHUTTER_STATE'][indices]\n", " srcx = hdul['SHUTTER_INFO'].data['ESTIMATED_SOURCE_IN_SHUTTER_X'][indices]\n", " srcy = hdul['SHUTTER_INFO'].data['ESTIMATED_SOURCE_IN_SHUTTER_Y'][indices]\n", " dith = hdul['SHUTTER_INFO'].data['DITHER_POINT_INDEX'][indices]\n", " primary = hdul['SHUTTER_INFO'].data['PRIMARY_SOURCE'][indices]\n", "\n", " print(slitlets)\n", " print(bkgd)\n", " print(primary)\n", "\n", " # construct a new table with just these values\n", " tabcol1 = fits.Column(name='SLITLET_ID', format='I', array=slitlets)\n", " tabcol2 = fits.Column(name='MSA_METADATA_ID', format='I', array=metaids)\n", " tabcol3 = fits.Column(name='SHUTTER_QUADRANT', format='I', array=quads)\n", " tabcol4 = fits.Column(name='SHUTTER_ROW', format='I', array=rows)\n", " tabcol5 = fits.Column(name='SHUTTER_COLUMN', format='I', array=cols)\n", " tabcol6 = fits.Column(name='SOURCE_ID', format='I', array=sources)\n", " tabcol7 = fits.Column(name='BACKGROUND', format='A', array=bkgd)\n", " tabcol8 = fits.Column(name='SHUTTER_STATE', format='4A', array=state)\n", " tabcol9 = fits.Column(name='ESTIMATED_SOURCE_IN_SHUTTER_X', format='E', array=srcx)\n", " tabcol10 = fits.Column(name='ESTIMATED_SOURCE_IN_SHUTTER_Y', format='E', array=srcy)\n", " tabcol11 = fits.Column(name='DITHER_POINT_INDEX', format='I', array=dith)\n", " tabcol12 = fits.Column(name='PRIMARY_SOURCE', format='1A', array=primary)\n", " hdul['SHUTTER_INFO'] = fits.BinTableHDU.from_columns([tabcol1, tabcol2, tabcol3, tabcol4, tabcol5,\n", " tabcol6, tabcol7, tabcol8, tabcol9, tabcol10,\n", " tabcol11, tabcol12], name='SHUTTER_INFO')\n", "\n", " # save to a new metafile\n", " hdul.writeto(output_dir+'newmetafile_1shutter_msa.fits', overwrite=True)\n", " hdul.close()" ] }, { "cell_type": "code", "execution_count": null, "id": "298f81a2", "metadata": {}, "outputs": [], "source": [ "# take a look at the contents of this new metafile\n", "\n", "metafile = output_dir+'newmetafile_1shutter_msa.fits'\n", "hdul = fits.open(metafile)\n", "\n", "print(hdul['SHUTTER_INFO'].data)\n", "\n", "hdul.close()" ] }, { "cell_type": "markdown", "id": "4b17eef3", "metadata": {}, "source": [ "Copy the rate files for the single exposure we've selected (corresponding to the first nod position) so that the pipeline products don't get overwritten. Also, we need to change the metafile header keyword to use the new file we created in the previous cell." ] }, { "cell_type": "code", "execution_count": null, "id": "fb2a0703", "metadata": {}, "outputs": [], "source": [ "if runflag:\n", " copyfile(output_dir+'nirspec_mossim_d1_g235m_nrs1_rate.fits',\n", " output_dir+'nirspec_mossim_d1_g235m_nrs1_1shutter_rate.fits')\n", " hdul = fits.open(output_dir+'nirspec_mossim_d1_g235m_nrs1_1shutter_rate.fits', mode='update')\n", " print(hdul[0].header['MSAMETFL'])\n", " hdul[0].header['MSAMETFL'] = 'newmetafile_1shutter_msa.fits'\n", " print(hdul[0].header['MSAMETFL'])\n", " hdul.close()\n", "\n", " copyfile(output_dir+'nirspec_mossim_d1_g235m_nrs2_rate.fits',\n", " output_dir+'nirspec_mossim_d1_g235m_nrs2_1shutter_rate.fits')\n", " hdul = fits.open(output_dir+'nirspec_mossim_d1_g235m_nrs2_1shutter_rate.fits', mode='update')\n", " print(hdul[0].header['MSAMETFL'])\n", " hdul[0].header['MSAMETFL'] = 'newmetafile_1shutter_msa.fits'\n", " print(hdul[0].header['MSAMETFL'])\n", " hdul.close()" ] }, { "cell_type": "code", "execution_count": null, "id": "bbb66194", "metadata": {}, "outputs": [], "source": [ "# reprocess through spec2\n", "\n", "if runflag:\n", " input_file = output_dir+'nirspec_mossim_d1_g235m_nrs1_1shutter_rate.fits'\n", " spec2 = Spec2Pipeline()\n", " spec2.save_results = True\n", " spec2.output_dir = output_dir\n", " spec2.master_background.save_background = True # output the master background spectrum\n", " result = spec2(input_file)\n", " \n", " input_file = output_dir+'nirspec_mossim_d1_g235m_nrs2_1shutter_rate.fits'\n", " spec2 = Spec2Pipeline()\n", " spec2.save_results = True\n", " spec2.output_dir = output_dir\n", " spec2.master_background.save_background = True # output the master background spectrum\n", " result = spec2(input_file)" ] }, { "cell_type": "markdown", "id": "07e641b9", "metadata": {}, "source": [ "Next, we will repeat the above, but now adding in one-shutter slitlets marked as being open to background. The pipeline will read this metafile and automatically identify that a master background needs to be constructed and subtracted from each slitlet. This simulation was not specifically designed with dedicated background slitlets, so we will create them by hand using background shutters in slitlets where the source flux was not detectable." ] }, { "cell_type": "code", "execution_count": null, "id": "fcbdf2d3", "metadata": {}, "outputs": [], "source": [ "if runflag:\n", "\n", " metafile = output_dir+'pointing_summary_n_metafile_msa.fits'\n", " hdul = fits.open(metafile)\n", " \n", " # pick the same source shutter as before\n", " # add shutters from other slitlets that are not contaminated by source flux to use as dedicated background slitlets\n", " indices = [255, 261, 264, 267, 288, 291, 294, 324, 327, 330, 360, 363, 366]\n", "\n", " slitlets = hdul['SHUTTER_INFO'].data['SLITLET_ID'][indices]\n", " metaids = hdul['SHUTTER_INFO'].data['MSA_METADATA_ID'][indices]\n", " quads = hdul['SHUTTER_INFO'].data['SHUTTER_QUADRANT'][indices]\n", " rows = hdul['SHUTTER_INFO'].data['SHUTTER_ROW'][indices]\n", " cols = hdul['SHUTTER_INFO'].data['SHUTTER_COLUMN'][indices]\n", " sources = hdul['SHUTTER_INFO'].data['SOURCE_ID'][indices]\n", " bkgd = hdul['SHUTTER_INFO'].data['BACKGROUND'][indices]\n", " state = hdul['SHUTTER_INFO'].data['SHUTTER_STATE'][indices]\n", " srcx = hdul['SHUTTER_INFO'].data['ESTIMATED_SOURCE_IN_SHUTTER_X'][indices]\n", " srcy = hdul['SHUTTER_INFO'].data['ESTIMATED_SOURCE_IN_SHUTTER_Y'][indices]\n", " dith = hdul['SHUTTER_INFO'].data['DITHER_POINT_INDEX'][indices]\n", " primary = hdul['SHUTTER_INFO'].data['PRIMARY_SOURCE'][indices]\n", "\n", " # set BACKGROUND=Y and SOURCE=N for all shutters we want to use to construct the master background\n", " bkgd[1:12] = 'Y'\n", " primary[1:12] = 'N'\n", " # give these different slitlet numbers so the pipeline will extract them separately\n", " slitlets = [31, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111]\n", "\n", " print(slitlets)\n", " print(bkgd)\n", " print(primary)\n", " \n", " # construct a new table with just these elements\n", " tabcol1 = fits.Column(name='SLITLET_ID', format='I', array=slitlets)\n", " tabcol2 = fits.Column(name='MSA_METADATA_ID', format='I', array=metaids)\n", " tabcol3 = fits.Column(name='SHUTTER_QUADRANT', format='I', array=quads)\n", " tabcol4 = fits.Column(name='SHUTTER_ROW', format='I', array=rows)\n", " tabcol5 = fits.Column(name='SHUTTER_COLUMN', format='I', array=cols)\n", " tabcol6 = fits.Column(name='SOURCE_ID', format='I', array=sources)\n", " tabcol7 = fits.Column(name='BACKGROUND', format='A', array=bkgd)\n", " tabcol8 = fits.Column(name='SHUTTER_STATE', format='4A', array=state)\n", " tabcol9 = fits.Column(name='ESTIMATED_SOURCE_IN_SHUTTER_X', format='E', array=srcx)\n", " tabcol10 = fits.Column(name='ESTIMATED_SOURCE_IN_SHUTTER_Y', format='E', array=srcy)\n", " tabcol11 = fits.Column(name='DITHER_POINT_INDEX', format='I', array=dith)\n", " tabcol12 = fits.Column(name='PRIMARY_SOURCE', format='1A', array=primary)\n", " hdul['SHUTTER_INFO'] = fits.BinTableHDU.from_columns([tabcol1, tabcol2, tabcol3, tabcol4, tabcol5,\n", " tabcol6, tabcol7, tabcol8, tabcol9, tabcol10,\n", " tabcol11, tabcol12], name='SHUTTER_INFO')\n", "\n", " # save to a new metafile\n", " hdul.writeto(output_dir+'newmetafile_1shutter+background_msa.fits', overwrite=True)\n", " hdul.close()" ] }, { "cell_type": "code", "execution_count": null, "id": "915f8641", "metadata": {}, "outputs": [], "source": [ "# take a look at the contents of this new metafile\n", "\n", "metafile = output_dir+'newmetafile_1shutter+background_msa.fits'\n", "hdul = fits.open(metafile)\n", "\n", "print(hdul['SHUTTER_INFO'].data)\n", "\n", "hdul.close()" ] }, { "cell_type": "markdown", "id": "5855f395", "metadata": {}, "source": [ "Once again, copy the rate files for the single exposure we've selected (corresponding to the first nod position) so that the pipeline products don't get overwritten. Also, we need to change the metafile header keyword to use the new file we created in the previous cell." ] }, { "cell_type": "code", "execution_count": null, "id": "c6f6c7c9", "metadata": {}, "outputs": [], "source": [ "if runflag:\n", " \n", " copyfile(output_dir+'nirspec_mossim_d1_g235m_nrs1_rate.fits',\n", " output_dir+'nirspec_mossim_d1_g235m_nrs1_1shutter_MBS_rate.fits')\n", " hdul = fits.open(output_dir+'nirspec_mossim_d1_g235m_nrs1_1shutter_MBS_rate.fits', mode='update')\n", " print(hdul[0].header['MSAMETFL'])\n", " hdul[0].header['MSAMETFL'] = 'newmetafile_1shutter+background_msa.fits'\n", " print(hdul[0].header['MSAMETFL'])\n", " hdul.close()\n", "\n", " copyfile(output_dir+'nirspec_mossim_d1_g235m_nrs2_rate.fits',\n", " output_dir+'nirspec_mossim_d1_g235m_nrs2_1shutter_MBS_rate.fits')\n", " hdul = fits.open(output_dir+'nirspec_mossim_d1_g235m_nrs2_1shutter_MBS_rate.fits', mode='update')\n", " print(hdul[0].header['MSAMETFL'])\n", " hdul[0].header['MSAMETFL'] = 'newmetafile_1shutter+background_msa.fits'\n", " print(hdul[0].header['MSAMETFL'])\n", " hdul.close()" ] }, { "cell_type": "markdown", "id": "604a0459", "metadata": {}, "source": [ "Now we will rerun the calwebb_spec2 pipeline using the new metafile. The processing proceeds as follows:\n", "\n", "* All background slitlets are processed first, calibrated as an extended source. If there are multiple background slitlets, all the extracted 1D spectra are averaged together to create a single master background spectrum.\n", "* Each source slitlet is processed through part of the spec2 pipeline up to the master_background step (including the assign_wcs, imprint, msaflagopen, extract_2d, and srctype steps).\n", "* The 1D master background is then resampled into the unrectified 2D space of the source slitlet, and the various throughput-related corrections (e.g., flat_field, barshadow, etc) are *removed*.\n", "* The 2D background is subtracted from the 2D source spectrum.\n", "* The remaining spec2 steps are run on the subtracted source spectrum." ] }, { "cell_type": "code", "execution_count": null, "id": "978e0a31", "metadata": {}, "outputs": [], "source": [ "if runflag:\n", " input_file = output_dir+'nirspec_mossim_d1_g235m_nrs1_1shutter_MBS_rate.fits'\n", " spec2 = Spec2Pipeline()\n", " spec2.save_results = True\n", " spec2.output_dir = output_dir\n", " spec2.master_background.save_background = True # output the master background spectrum\n", " result = spec2(input_file)\n", " \n", " input_file = output_dir+'nirspec_mossim_d1_g235m_nrs2_1shutter_MBS_rate.fits'\n", " spec2 = Spec2Pipeline()\n", " spec2.save_results = True\n", " spec2.output_dir = output_dir\n", " spec2.master_background.save_background = True # output the master background spectrum\n", " result = spec2(input_file)" ] }, { "cell_type": "code", "execution_count": null, "id": "3ef5054f", "metadata": {}, "outputs": [], "source": [ "# show the results\n", "\n", "calfile = output_dir+'nirspec_mossim_d1_g235m_nrs1_1shutter_MBS_cal.fits'\n", "cal = datamodels.open(calfile)\n", "root = calfile[:-9]\n", "s2d = datamodels.open(root+'_s2d.fits') # this contains the calibrated *rectified* 2D spectrum\n", "x1d = datamodels.open(root+'_x1d.fits') # this contains the aperture-extracted 1D spectrum\n", "mb2d = datamodels.open(root+'_masterbg2d.fits') # this contains the 2D master background spectrum\n", "mbx1d = datamodels.open(root+'_masterbg1d.fits') # this contains the 1D master background spectrum\n", "\n", "# first, display the master background spectrum that was created\n", "# the combined 1D spectrum\n", "mbx1dwave = mbx1d.spec[0].spec_table.WAVELENGTH\n", "mbx1dflux = mbx1d.spec[0].spec_table.FLUX\n", "fig = plt.figure(figsize=(19, 8))\n", "plt.plot(mbx1dwave, mbx1dflux)\n", "plt.xlabel('wavelength (microns)')\n", "plt.ylabel('flux (Jy)')\n", "plt.title('1D master background')\n", "plt.ylim(-2.e3, 1.5e4)\n", "plt.show()\n", "\n", "# 2D MB (resampled into the unrectified 2D space of slitlet 31)\n", "mb2dsci = mb2d.slits[0].data\n", "show_image(mb2dsci, 0.0, 0.02, aspect=10., scale='linear', units='MJy')\n", "\n", "# now show the source spectrum\n", "for i, slit in enumerate(cal.slits):\n", " # let's look at the subtracted source\n", " if (slit.name == '31'):\n", " print(slit.name)\n", " \n", " calsci = slit.data # contains the pixel data from the cal file (SCI extension)\n", " s2dsci = s2d.slits[i].data # contains the pixel data from the s2d file\n", " \n", " # determine the wavelength scale of the s2d data for plotting purposes\n", " # get the data model WCS object\n", " wcsobj = s2d.slits[i].meta.wcs\n", " y, x = np.mgrid[:s2dsci.shape[0], : s2dsci.shape[1]] # grid of pixel x,y indices\n", " det2sky = wcsobj.get_transform('detector', 'world') # the coordinate transform from detector space (pixels) to sky (RA, DEC in degrees)\n", " ra, dec, s2dwave = det2sky(x, y) # RA, Dec, wavelength (microns) for each pixel\n", " s2dwaves = s2dwave[0, :] # only need a single row of values since this is the rectified spectrum\n", " xtint = np.arange(100, s2dsci.shape[1], 100)\n", " xtlab = np.round(s2dwaves[xtint], 2) # wavelength labels for the x-axis\n", " \n", " # get wavelength & flux from the x1d data model\n", " x1dwave = x1d.spec[i].spec_table.WAVELENGTH\n", " x1dflux = x1d.spec[i].spec_table.FLUX\n", " \n", " # plot the unrectified calibrated 2D spectrum\n", " show_image(calsci, -0.01, 0.01, aspect=10., scale='linear', units='MJy')\n", " \n", " # plot the rectified 2D spectrum\n", " show_image(s2dsci, -0.01, 0.01, aspect=10., scale='linear', units='MJy')\n", " plt.xticks(xtint, xtlab)\n", " plt.xlabel('wavelength (microns)')\n", " \n", " # plot the 1D extracted spectrum\n", " fig = plt.figure(figsize=(19, 8))\n", " plt.plot(x1dwave, x1dflux)\n", " plt.xlabel('wavelength (microns)')\n", " plt.ylabel('flux (Jy)')\n", " plt.title('background-subtracted source')\n", " plt.ylim(-1.e4, 1.e6)\n", " plt.show()" ] }, { "cell_type": "markdown", "id": "9cb83547", "metadata": {}, "source": [ "Finally, let's compare all three x1d products generated for this source: the unsubtracted version, the nod-subtracted product, and the master background-subtracted product." ] }, { "cell_type": "code", "execution_count": null, "id": "76bc96e2", "metadata": {}, "outputs": [], "source": [ "x1dnosub = datamodels.open(output_dir+'nirspec_mossim_d1_g235m_nrs1_1shutter_x1d.fits')\n", "x1dnodsub = datamodels.open(output_dir+'nirspec_mossim_d1_g235m_nrs1_x1d.fits')\n", "x1dmbsub = datamodels.open(output_dir+'nirspec_mossim_d1_g235m_nrs1_1shutter_MBS_x1d.fits')\n", "\n", "x1dnosubwave = x1dnosub.spec[0].spec_table.WAVELENGTH\n", "x1dnosubflux = x1dnosub.spec[0].spec_table.FLUX\n", "\n", "# need to loop through all the spectra for the nod-subtracted case to find source 2140\n", "for spec in x1dnodsub.spec:\n", " if (spec.source_id == 2140):\n", " x1dnodsubwave = spec.spec_table.WAVELENGTH\n", " x1dnodsubflux = spec.spec_table.FLUX\n", "\n", "x1dmbsubwave = x1dmbsub.spec[12].spec_table.WAVELENGTH\n", "x1dmbsubflux = x1dmbsub.spec[12].spec_table.FLUX\n", " \n", "# plot the 1D spectra\n", "fig = plt.figure(figsize=(19, 8))\n", "plt.plot(x1dnosubwave, x1dnosubflux)\n", "plt.plot(x1dnodsubwave, x1dnodsubflux)\n", "plt.title('results of nod subtraction')\n", "plt.ylim(0., 6e4)\n", "plt.xlabel('wavelength (microns)')\n", "plt.ylabel('flux (Jy)')\n", "plt.show()\n", "\n", "# plot the 1D extracted spectrum\n", "fig = plt.figure(figsize=(19, 8))\n", "plt.plot(x1dnosubwave, x1dnosubflux)\n", "plt.plot(x1dmbsubwave, x1dmbsubflux)\n", "plt.title('results of master background subtraction')\n", "plt.ylim(0., 6e4)\n", "plt.xlabel('wavelength (microns)')\n", "plt.ylabel('flux (Jy)')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "95cfcac8", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.12" } }, "nbformat": 4, "nbformat_minor": 5 }