diff --git a/_toc.yml b/_toc.yml index b6a95e131..6ed804370 100644 --- a/_toc.yml +++ b/_toc.yml @@ -53,10 +53,6 @@ parts: title: Improving Accuracy of WCS of Pure Parallel Dataset - caption: MIRI chapters: - - file: notebooks/MIRI/MRS_Mstar_analysis/JWST_Mstar_dataAnalysis_runpipeline.ipynb - title: MRS Mstar - Run Pipeline - - file: notebooks/MIRI/MRS_Mstar_analysis/JWST_Mstar_dataAnalysis_PointSourceDetectorBasedExtraction.ipynb - title: MRS Mstar - Optimal Extraction - file: notebooks/MIRI/MRS_Mstar_analysis/JWST_Mstar_dataAnalysis_analysis.ipynb title: MRS Mstar - Data Analysis - file: notebooks/MIRI/MIRI_IFU_YSOs_in_the_LMC/isha_nayak_ysos_in_the_lmc.ipynb diff --git a/notebooks/MIRI/MRS_Mstar_analysis/JWST_Mstar_dataAnalysis_PointSourceDetectorBasedExtraction.ipynb b/notebooks/MIRI/MRS_Mstar_analysis/JWST_Mstar_dataAnalysis_PointSourceDetectorBasedExtraction.ipynb deleted file mode 100644 index 4858e432d..000000000 --- a/notebooks/MIRI/MRS_Mstar_analysis/JWST_Mstar_dataAnalysis_PointSourceDetectorBasedExtraction.ipynb +++ /dev/null @@ -1,1368 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": { - "slideshow": { - "slide_type": "slide" - } - }, - "source": [ - "# MRS Detector Based Extraction" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "**Use case:** This notebook illustrates a spectral extraction technique from detector image plane using a PSF model and the pixel variance. This extraction is considered more advanced than what is currently implemented in the JWST pipeline and is expected to be fully integrated into the pipeline in the future.
\n", - "**Data:** Simulated [MIRI MRS](https://jwst-docs.stsci.edu/mid-infrared-instrument/miri-observing-modes/miri-medium-resolution-spectroscopy) spectrum of AGB star.
\n", - "**Tools:** jwst, astropy, numpy, scipy.
\n", - "**Cross-intrument:** No.
\n", - "**Documentation:** This notebook is part of a STScI's larger [post-pipeline Data Analysis Tools Ecosystem](https://jwst-docs.stsci.edu/jwst-post-pipeline-data-analysis) and can be [downloaded](https://github.com/spacetelescope/dat_pyinthesky/tree/main/jdat_notebooks/MRS_Mstar_analysis) directly from the [JDAT Notebook Github directory](https://github.com/spacetelescope/jdat_notebooks).
\n", - "**Source of Simulations:** [MIRISim](https://www.stsci.edu/jwst/science-planning/proposal-planning-toolbox/mirisim)
\n", - "**Pipeline Version:** [JWST Pipeline](https://jwst-docs.stsci.edu/jwst-data-reduction-pipeline)
\n", - "\n", - "**Note**: This notebook includes MIRI simulated data cubes obtained using MIRISim (https://wiki.miricle.org//bin/view/Public/MIRISim_Public) and run through the JWST pipeline (https://jwst-pipeline.readthedocs.io/en/latest/) of point sources with spectra representative of late M type stars.\n", - "\n", - "## Introduction\n", - "\n", - "This notebook analyzes one star represented by a dusty SED corresponding to the ISO SWS spectrum of W Per from Kraemer et al. (2002) and Sloan et al. (2003) to cover the MRS spectral range 5-28 microns. Analysis of JWST spectral cubes requires extracting spatial-spectral features of interest and measuring their attributes.\n", - "\n", - "This is the first notebook, which will process the data and automatically detect and extract spectra for the point source. The workflow automatically detects sources on the detector image plane to extract the spectrum of the point source. Then it will read in the spectra generated at Stage 3 of the JWST pipeline.\n", - "\n", - "## Motivation to the Detector Based Extraction Technique\n", - "\n", - "At present, one-dimensional spectra produced by the JWST pipeline are\n", - " constructed exclusively from the combined and rectified 3-dimensional\n", - "data cubes. These data cubes are treated as a collection of images at\n", - "various wavelength planes, and simple aperture-based (or effective\n", - "PSF-weighted) techniques are used to recover the source intensity as a\n", - "function of wavelength in a manner akin to standard imaging PSF\n", - "photometry. However, this method of extracting point source spectra\n", - "from data cubes is not ideal; fundamentally, a source with no spatial\n", - "structure should not undergo spatial reprocessing prior to extraction of\n", - "spectral information. Indeed, large-scale IFU surveys targeting point\n", - "sources rather than extended sources (e.g., SDSS IV/MaStar, R. Yan et\n", - "al. 2019) have relied on detector-based spectral extraction rather than\n", - "the construction of 3d data cubes.\n", - " \n", - "For the MIRI MRS in particular, such detector-based extraction offers\n", - "the opportunity to correct for the following effects which would be\n", - "difficult or impossible to handle once the observational data has been\n", - "rectified into 3-dimensional cubes:\n", - "\n", - "1) Given the significant undersampling of the MRS, offsets in spaxel\n", - "sampling locations as a function of wavelength can alias spatial\n", - "structure from the point spread function into spectral artifacts that\n", - "bias extracted 1D spectra. 3D cube building is designed in a manner to\n", - "mitigate these artifacts, but the act of resampling ensures that 1D\n", - "spectra of unresolved sources extracted from 3D cubes can never be as\n", - "reliable as 1D spectra extracted directly from detector-level data\n", - "without the intermediate resampling step.\n", - "\n", - "2) Resampled pixels in the data cubes have significant covariance with\n", - "adjacent/nearby pixels. Treating this covariance properly in the\n", - "construction of extracted 1D spectra with flux and variance information\n", - "is challenging, and not handled at present by the cube-based spectral\n", - "extraction.\n", - "\n", - "3) Pipeline calibrations of the instrument throughput and wavelength\n", - "calibration are based on sources that fill a given IFU slice, whereas\n", - "point sources produce a non-uniform response depending upon the location\n", - "of the source centroid within a given IFU slice. For instance, given\n", - "the degeneracy on the detector between the spatial across-slice and\n", - "spectral dispersion directions, the wavelength solution for point\n", - "sources will differ slightly than for extended sources in a manner that\n", - "depends on the source location within the slice. This effect cannot be\n", - "recovered once a 3D cube has been constructed that combines data from\n", - "many slices (and indeed, multiple dithered exposures).\n", - "\n", - "4) The MRS has a spectral leak in which 2nd order light from 6 microns\n", - "(Ch1B) is visible at a few percent transmission in the 12 micron\n", - "observations (Ch3A). This effect is not possible to correct in general,\n", - "since the Ch3A field of view is larger than that of Ch1B. However, for\n", - "the specific case of a point source within the common field of view of\n", - "all channels it is possible to use the Ch1B spectrum to predict and\n", - "subtract the second-order signal from the Ch3A data.\n", - "\n", - "## Implementation of the Algorithm\n", - "\n", - "In this notebook, we implement a method for achieving MIRI MRS\n", - "detector-based point source extraction. The formalism behind the method\n", - "follows the optimal spectral extraction formalism by Keith Horne 1986.\n", - "Extracting a 1d integrated spectrum directly from the MRS detectors\n", - "depends on (a) a detailed knowledge of the detector PSF, and (b) the\n", - "pixel variance. The MRS detector PSF differs considerably from the\n", - "telescope PSF due to the scattering of the photons inside the detector\n", - "substrate. How well the detector PSF is understood/modeled/measured is\n", - "key in this instance and should not be understated. For the pixel\n", - "variance, the slope-fitting uncertainty from the pixel integration\n", - "ramp(s) is used. The combination of the detector PSF and pixel variance\n", - "is used as a weighing function to weigh the signal in each detector pixel.\n", - "\n", - "In its most basic form, the notebook code assumes that there is only one\n", - "point source in the MRS field of view. Using a Gaussian-fitting routine\n", - "at detector-level, the centroid of the point source is determined with a\n", - "sub-pixel accuracy. The centroid is used to project the known detector\n", - "PSF model directly onto the point source signal. The centroid is also\n", - "used to derive the throughput and wavelength correction due to the\n", - "location of the point source with respect to the imaging slicer of the\n", - "MRS. The throughput and wavelength correction are applied before the\n", - "spectral extraction step.\n", - "\n", - "The spectral extraction step integrates the point source flux in a\n", - "pre-defined list of wavelength bins. All detector pixels contributing to\n", - "a wavelength bin are weighed by the projected PSF model -and- the pixel\n", - "variance. In case the extracted spectral band is MRS band 3A (around 12\n", - "microns) then the following steps are taken:\n", - "- The band 1B (~6 micron) detector image is multiplied by a spectral\n", - "leak transmission curve.\n", - "- A 1d “spectral leak” integrated point source spectrum is extracted\n", - "from band 1B.\n", - "- The spectral leak spectrum is subtracted from the band 3A integrated\n", - "point source spectrum.\n", - "\n", - "Future improvements to the code will allow users to input MRS point\n", - "source centroid locations to address the case of crowded fields, barely\n", - "resolved blended point sources, and point sources located within a\n", - "semi-extended environment with spatial structure. In addition, the fine\n", - "centroiding, currently performed using a Gaussian distribution, can be\n", - "improved by striving for residual minimisation with a master\n", - "detector-based PSF model. Such a master PSF model will be derived and\n", - "refined throughout commissioning and Cycle 1 calibration. \n", - "\n" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "slideshow": { - "slide_type": "slide" - } - }, - "source": [ - "## Import Packages" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "slideshow": { - "slide_type": "fragment" - } - }, - "outputs": [], - "source": [ - "# Import useful python packages\n", - "import numpy as np\n", - "from scipy.interpolate import interp1d\n", - "\n", - "# Import packages to display images in the notebook\n", - "from matplotlib.lines import Line2D\n", - "import matplotlib.pyplot as plt \n", - "%matplotlib inline\n", - "\n", - "# Set general plotting options\n", - "params = {'legend.fontsize': '18', 'axes.labelsize': '18', \n", - " 'axes.titlesize': '18', 'xtick.labelsize': '18', \n", - " 'ytick.labelsize': '18', 'lines.linewidth': 2, 'axes.linewidth': 2, 'animation.html': 'html5'}\n", - "plt.rcParams.update(params)\n", - "plt.rcParams.update({'figure.max_open_warning': 0})\n", - "\n", - "from IPython.display import display, HTML\n", - "display(HTML(\"\"))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import warnings\n", - "warnings.simplefilter('ignore')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "jupyter": { - "outputs_hidden": false - }, - "pycharm": { - "name": "#%%\n" - } - }, - "outputs": [], - "source": [ - "# Import astropy packages \n", - "from astropy import units as u\n", - "from astropy.io import ascii\n", - "from astropy.wcs import WCS\n", - "from astropy.table import Table, vstack\n", - "from astropy.stats import sigma_clipped_stats\n", - "from astropy.nddata import StdDevUncertainty\n", - "from astropy.io import fits # added by BAS on 8 April 2021\n", - "from astropy.utils.data import get_pkg_data_filename\n", - "from astropy.utils.data import download_file\n", - "\n", - "# To find stars in the MRS spectralcubes and do aperture photometry\n", - "from photutils import DAOStarFinder, CircularAperture\n", - "\n", - "# To deal with 1D spectrum\n", - "from specutils import Spectrum1D\n", - "from specutils.fitting import fit_generic_continuum\n", - "from specutils.manipulation import box_smooth, extract_region, SplineInterpolatedResampler\n", - "from specutils.analysis import line_flux, centroid, equivalent_width\n", - "from specutils.spectra import SpectralRegion\n", - "from specutils import SpectrumList\n", - "\n", - "# To make nice plots with WCS axis\n", - "# import aplpy\n", - "\n", - "# To fit a curve to the data\n", - "from scipy.optimize import curve_fit\n", - "\n", - "# To download data\n", - "import tarfile\n", - "import urllib.request\n", - "from pathlib import Path" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# import necessary functions from PointSourceDetectorBasedExtractionFuncs and create_distortionMaps\n", - "from PointSourceDetectorBasedExtractionFuncs import mrs_aux, point_source_centroiding, evaluate_psf_cdp\n", - "\n", - "# evaluate distortion maps on the 2D detector image plane\n", - "from create_distortionMaps import d2cMapping" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Set paths to the Data and Outputs\n", - "\n", - "Use MIRISim JWST pipeline processed data in future iterations." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import os\n", - "os.environ['CRDS_PATH'] = os.environ['HOME']+'/crds_cache'\n", - "os.environ['CRDS_SERVER_URL'] = 'https://jwst-crds-pub.stsci.edu'" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Set JWST pipeline paths\n", - "import os\n", - "os.environ['CRDS_PATH'] = str(Path(os.environ['HOME']) / 'crds_cache')\n", - "os.environ['CRDS_SERVER_URL'] = 'https://jwst-crds-pub.stsci.edu'\n", - "\n", - "# import JWST pipeline\n", - "from jdaviz.app import Application\n", - "from jwst.pipeline import Detector1Pipeline\n", - "from jwst.pipeline import Spec2Pipeline\n", - "from jwst.pipeline import Spec3Pipeline\n", - "from jwst.extract_1d import Extract1dStep\n", - "from jwst.associations.lib.rules_level3_base import DMS_Level3_Base\n", - "from jwst.associations import asn_from_list\n", - "from photutils import aperture_photometry\n", - "\n", - "import asdf\n", - "import crds\n", - "import json\n", - "import glob" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tar_location = Path('./')\n", - "tar_files = [\n", - " '20210413_120546_mirisim.tar.gz', \n", - " '20210413_123047_mirisim.tar.gz',\n", - " '20210413_125354_mirisim.tar.gz'\n", - "]\n", - "tar_paths_full = [tar_location / f for f in tar_files]\n", - "\n", - "tar_urls = [\n", - " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/20210413_120546_mirisim.tar.gz',\n", - " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/20210413_123047_mirisim.tar.gz',\n", - " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/20210413_125354_mirisim.tar.gz'\n", - "]\n", - "\n", - "\n", - "for p, url in zip(tar_paths_full, tar_urls):\n", - " # check for un-tarred version of file\n", - " tmp = p.parent / p.stem\n", - " if os.path.exists(tmp.stem):\n", - " print(f\"{p.name} Directory Exists.\")\n", - " else:\n", - " # check for tarred version of file and un-tar\n", - " if not os.path.exists(p):\n", - " # download file if absent\n", - " print(f\"Downloading {p.name}...\") \n", - " urllib.request.urlretrieve(url, p)\n", - " \n", - " tar = tarfile.open(p, \"r:gz\")\n", - " tar.extractall()\n", - " tar.close()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Create Directories If They Don't Exist\n", - "if os.path.isdir(\"./miri_devel/cdp_data/DISTORTION\"):\n", - " print(\"MIRI Reference File Directory Exists.\")\n", - "else:\n", - " os.makedirs(\"./miri_devel/cdp_data/DISTORTION\")\n", - " print(\"Made Directory ./miri_devel/cdp_data/DISTORTION\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "data_dir = Path('./')\n", - "data_paths = [\n", - " 'spectral_leak_transmission.npy',\n", - " 'W_Per_spectrum_sws.fit',\n", - " 'jwst_miri_photom_0064.fits',\n", - " 'jwst_miri_photom_0060.fits',\n", - " 'jwst_miri_photom_0052.fits',\n", - " 'miri_devel/cdp_data/MIRI_FM_MIRIFUSHORT_2SHORT_PSF_07.02.00.fits',\n", - " 'miri_devel/cdp_data/MIRI_FM_MIRIFUSHORT_2MEDIUM_PSF_07.02.00.fits',\n", - " 'miri_devel/cdp_data/MIRI_FM_MIRIFUSHORT_2LONG_PSF_07.02.00.fits',\n", - " 'miri_devel/cdp_data/MIRI_FM_MIRIFUSHORT_1SHORT_PSF_07.02.00.fits',\n", - " 'miri_devel/cdp_data/MIRI_FM_MIRIFUSHORT_1MEDIUM_PSF_07.02.00.fits',\n", - " 'miri_devel/cdp_data/MIRI_FM_MIRIFUSHORT_1LONG_PSF_07.02.00.fits',\n", - " 'miri_devel/cdp_data/MIRI_FM_MIRIFULONG_3SHORT_PSF_07.02.00.fits',\n", - " 'miri_devel/cdp_data/MIRI_FM_MIRIFULONG_3MEDIUM_PSF_07.02.00.fits',\n", - " 'miri_devel/cdp_data/MIRI_FM_MIRIFULONG_3LONG_PSF_07.02.00.fits',\n", - " 'miri_devel/cdp_data/MIRI_FM_MIRIFULONG_4SHORT_PSF_07.02.00.fits',\n", - " 'miri_devel/cdp_data/MIRI_FM_MIRIFULONG_4MEDIUM_PSF_07.02.00.fits',\n", - " 'miri_devel/cdp_data/MIRI_FM_MIRIFULONG_4LONG_PSF_07.02.00.fits',\n", - " 'miri_devel/cdp_data/DISTORTION/MIRI_FM_MIRIFULONG_34LONG_DISTORTION_8B.05.03.fits',\n", - " 'miri_devel/cdp_data/DISTORTION/MIRI_FM_MIRIFULONG_34LONG_DISTORTION_8B.05.02.fits',\n", - " 'miri_devel/cdp_data/DISTORTION/MIRI_FM_MIRIFULONG_34MEDIUM_DISTORTION_8B.05.02.fits',\n", - " 'miri_devel/cdp_data/DISTORTION/MIRI_FM_MIRIFULONG_34MEDIUM_DISTORTION_8B.05.03.fits',\n", - " 'miri_devel/cdp_data/DISTORTION/MIRI_FM_MIRIFULONG_34SHORT_DISTORTION_8B.05.02.fits',\n", - " 'miri_devel/cdp_data/DISTORTION/MIRI_FM_MIRIFULONG_34SHORT_DISTORTION_8B.05.03.fits',\n", - " 'miri_devel/cdp_data/DISTORTION/MIRI_FM_MIRIFUSHORT_12LONG_DISTORTION_8B.05.02.fits',\n", - " 'miri_devel/cdp_data/DISTORTION/MIRI_FM_MIRIFUSHORT_12MEDIUM_DISTORTION_8B.05.02.fits',\n", - " 'miri_devel/cdp_data/DISTORTION/MIRI_FM_MIRIFUSHORT_12SHORT_DISTORTION_8B.05.02.fits'\n", - "]\n", - "\n", - "data_paths_full = [data_dir / p for p in data_paths]\n", - "data_urls = [\n", - " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/spectral_leak_transmission.npy',\n", - " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/W_Per_spectrum_sws.fit',\n", - " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/jwst_miri_photom_0064.fits',\n", - " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/jwst_miri_photom_0060.fits',\n", - " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/jwst_miri_photom_0052.fits',\n", - " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/miri_devel/cdp_data/MIRI_FM_MIRIFUSHORT_2SHORT_PSF_07.02.00.fits',\n", - " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/miri_devel/cdp_data/MIRI_FM_MIRIFUSHORT_2MEDIUM_PSF_07.02.00.fits',\n", - " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/miri_devel/cdp_data/MIRI_FM_MIRIFUSHORT_2LONG_PSF_07.02.00.fits',\n", - " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/miri_devel/cdp_data/MIRI_FM_MIRIFUSHORT_1SHORT_PSF_07.02.00.fits',\n", - " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/miri_devel/cdp_data/MIRI_FM_MIRIFUSHORT_1MEDIUM_PSF_07.02.00.fits',\n", - " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/miri_devel/cdp_data/MIRI_FM_MIRIFUSHORT_1LONG_PSF_07.02.00.fits',\n", - " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/miri_devel/cdp_data/MIRI_FM_MIRIFULONG_3SHORT_PSF_07.02.00.fits',\n", - " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/miri_devel/cdp_data/MIRI_FM_MIRIFULONG_3MEDIUM_PSF_07.02.00.fits',\n", - " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/miri_devel/cdp_data/MIRI_FM_MIRIFULONG_3LONG_PSF_07.02.00.fits',\n", - " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/miri_devel/cdp_data/MIRI_FM_MIRIFULONG_4SHORT_PSF_07.02.00.fits',\n", - " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/miri_devel/cdp_data/MIRI_FM_MIRIFULONG_4MEDIUM_PSF_07.02.00.fits',\n", - " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/miri_devel/cdp_data/MIRI_FM_MIRIFULONG_4LONG_PSF_07.02.00.fits',\n", - " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/miri_devel/cdp_data/DISTORTION/MIRI_FM_MIRIFULONG_34LONG_DISTORTION_8B.05.03.fits',\n", - " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/miri_devel/cdp_data/DISTORTION/MIRI_FM_MIRIFULONG_34LONG_DISTORTION_8B.05.02.fits',\n", - " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/miri_devel/cdp_data/DISTORTION/MIRI_FM_MIRIFULONG_34MEDIUM_DISTORTION_8B.05.02.fits',\n", - " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/miri_devel/cdp_data/DISTORTION/MIRI_FM_MIRIFULONG_34MEDIUM_DISTORTION_8B.05.03.fits',\n", - " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/miri_devel/cdp_data/DISTORTION/MIRI_FM_MIRIFULONG_34SHORT_DISTORTION_8B.05.02.fits',\n", - " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/miri_devel/cdp_data/DISTORTION/MIRI_FM_MIRIFULONG_34SHORT_DISTORTION_8B.05.03.fits',\n", - " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/miri_devel/cdp_data/DISTORTION/MIRI_FM_MIRIFUSHORT_12LONG_DISTORTION_8B.05.02.fits',\n", - " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/miri_devel/cdp_data/DISTORTION/MIRI_FM_MIRIFUSHORT_12MEDIUM_DISTORTION_8B.05.02.fits',\n", - " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/miri_devel/cdp_data/DISTORTION/MIRI_FM_MIRIFUSHORT_12SHORT_DISTORTION_8B.05.02.fits'\n", - "] \n", - "\n", - "for p, url in zip(data_paths_full, data_urls):\n", - " if os.path.exists(p):\n", - " print(f\"{p.name} exists.\")\n", - " else:\n", - " print(f\"Downloading {p.name}...\")\n", - " urllib.request.urlretrieve(url, p)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Run Pipeline\n", - "\n", - "The various [stages of the pipeline](https://jwst-pipeline.readthedocs.io/en/latest/jwst/pipeline/main.html#pipelines) can be [run within Python](https://jwst-pipeline.readthedocs.io/en/latest/jwst/introduction.html#running-from-within-python). For a more in depth tutorial on running the pipelines, check out the [JWebbinars](https://www.stsci.edu/jwst/science-execution/jwebbinars)." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "# Execute calwebb_detector1 pipeline on raw simulation output. This will overwrite previous reductions.\n", - "\n", - "allshortfiles = glob.glob('20210413_*_mirisim/det_images/*MIRIFUSHORT*fits')\n", - "alllongfiles = glob.glob('20210413_*_mirisim/det_images/*MIRIFULONG*fits')\n", - " \n", - "pipe1short = Detector1Pipeline()\n", - "\n", - "# run calwebb_detector1 on the MIRIFUSHORT data separate from MIRIFULONG data, as it saves time this way\n", - "for shortfile in allshortfiles:\n", - " print(shortfile)\n", - " baseshort, remaindershort = shortfile.split('.')\n", - " \n", - " # If you run your own simulations, you will need to update these hardcoded files.\n", - " beforestuffshort, dateafterstuffshort = shortfile.split('20210413_') \n", - " datestringshort, afterstuffshort = dateafterstuffshort.split('_mirisim')\n", - " \n", - " pipe1short.refpix.skip = True\n", - " pipe1short.output_file = baseshort + datestringshort\n", - " \n", - " pipe1short.run(shortfile)\n", - "\n", - "pipe1long = Detector1Pipeline()\n", - "\n", - "for longfile in alllongfiles:\n", - " print(longfile)\n", - " baselong, remainderlong = longfile.split('.')\n", - " \n", - " # If you run your own simulations, you will need to update these hardcoded files.\n", - " beforestufflong, dateafterstufflong = longfile.split('20210413_')\n", - " datestringlong, afterstufflong = dateafterstufflong.split('_mirisim')\n", - " \n", - " pipe1long.refpix.skip = True\n", - " pipe1long.output_file = baselong + datestringlong\n", - " \n", - " pipe1long.run(longfile)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Execute calwebb_spec2 pipeline. This will overwrite previous reductions.\n", - "\n", - "# All the local calwebb_detector1 files\n", - "allshortfiles2 = glob.glob('det_image_*_MIRIFUSHORT_*_rate.fits')\n", - "alllongfiles2 = glob.glob('det_image_*_MIRIFULONG_*_rate.fits')\n", - "\n", - "for short2file in allshortfiles2:\n", - " print(short2file)\n", - " pipe2short = Spec2Pipeline()\n", - " base2short, remainder2short = short2file.split('.')\n", - " \n", - " pipe2short.straylight.skip = True\n", - " \n", - " # If you run your own simulations, you will need to update these hardcoded files.\n", - " if (short2file == 'det_image_seq1_MIRIFUSHORT_12LONGexp1125354_rate.fits'):\n", - " print('this one will have the level 2b cube built')\n", - " else:\n", - " pipe2short.cube_build.skip = True\n", - " pipe2short.extract_1d.skip = True\n", - " pipe2short.output_file = base2short\n", - " \n", - " pipe2short.run(short2file)\n", - "\n", - "for long2file in alllongfiles2:\n", - " print(long2file)\n", - " pipe2long = Spec2Pipeline()\n", - " base2long, remainder2long = long2file.split('.')\n", - " \n", - " pipe2long.straylight.skip = True\n", - " # If you run your own simulations, you will need to update these hardcoded files.\n", - " if (long2file == 'det_image_seq1_MIRIFULONG_34SHORTexp1120546_rate.fits'):\n", - " print('this one will have the level 2b cube built')\n", - " else:\n", - " pipe2long.cube_build.skip = True\n", - " pipe2long.extract_1d.skip = True\n", - " pipe2long.output_file = base2long\n", - " \n", - " pipe2long.run(long2file)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# If you run your own simulations, you will need to update these hardcoded files.\n", - "l_cube_file = 'det_image_seq1_MIRIFULONG_34SHORTexp1120546_s3d.fits'\n", - "s_cube_file = 'det_image_seq1_MIRIFUSHORT_12LONGexp1125354_s3d.fits'\n", - "\n", - "with fits.open(s_cube_file) as hdu_s_cube:\n", - " s_cube = hdu_s_cube['SCI'].data\n", - " s_med_cube = np.zeros((s_cube.shape[1], s_cube.shape[2]))\n", - " for a in range(s_cube.shape[1]):\n", - " for b in range(s_cube.shape[2]):\n", - " s_med_cube[a, b] = np.median(s_cube[:, a, b])\n", - "\n", - "mean, median, std = sigma_clipped_stats(s_med_cube, sigma=2.0)\n", - "\n", - "# Get a list of sources using a dedicated source detection algorithm\n", - "# Find sources at least 3* background (typically)\n", - "\n", - "daofind = DAOStarFinder(fwhm=2.0, threshold=3.*std)\n", - "sources = daofind(s_med_cube-median) \n", - "print(\"\\n Number of sources in field: \", len(sources))\n", - "\n", - "# Positions in pixels\n", - "positions = Table([sources['xcentroid'], sources['ycentroid']])\n", - "\n", - "# Convert to RA & Dec (ICRS)\n", - "peakpixval = np.zeros(len(sources['xcentroid']))\n", - "for count_s, _ in enumerate(sources):\n", - " peakpixval[count_s] = s_med_cube[int(np.round(sources['xcentroid'][count_s])), int(np.round(sources['ycentroid'][count_s]))]\n", - "print('peak pixel x =')\n", - "print(sources['xcentroid'][np.argmax(peakpixval)])\n", - "print('peak pixel y =')\n", - "print(sources['ycentroid'][np.argmax(peakpixval)])\n", - "\n", - "plt.imshow(s_med_cube, vmin=0, vmax=100)\n", - "plt.tight_layout()\n", - "plt.scatter(sources['xcentroid'], sources['ycentroid'], c=\"red\", marker=\"+\", s=50)\n", - "plt.scatter(sources['xcentroid'][np.argmax(peakpixval)], sources['ycentroid'][np.argmax(peakpixval)], c='black', marker='+', s=50)\n", - "plt.show()\n", - "\n", - "f0 = fits.open(s_cube_file)\n", - "w0 = WCS(f0[('sci', 1)].header, f0)\n", - "f0.close()\n", - "\n", - "radec = w0.all_pix2world([sources['xcentroid'][np.argmax(peakpixval)]], [sources['ycentroid'][np.argmax(peakpixval)]], [1], 1)\n", - "\n", - "# Take the brightest source flux and take that to be your primary point source for extraction\n", - "ra_ptsrc = radec[0][0]\n", - "dec_ptsrc = radec[1][0]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Due to the way the pipeline currently extracts Level3 data, you must update the headers to be centered on the point source of your choosing from the step above.\n", - "all_files = glob.glob('det_image_*_cal.fits')\n", - "targra = ra_ptsrc\n", - "targdec = dec_ptsrc\n", - "for thisfile in all_files:\n", - " base, remainder = thisfile.split('.')\n", - " outfilename = base + '_fix.' + remainder\n", - " print(outfilename)\n", - " \n", - " with fits.open(thisfile) as hduthis:\n", - " hduthis['SCI'].header['SRCTYPE'] = 'POINT'\n", - " hduthis[0].header['TARG_RA'] = targra \n", - " hduthis[0].header['TARG_DEC'] = targdec\n", - " hduthis.writeto(outfilename, overwrite=True)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# set up needed reference file(s) for spec3\n", - "\n", - "file_all_list = glob.glob('det_image_*_cal.fits')\n", - "\n", - "asnall = asn_from_list.asn_from_list(file_all_list, rule=DMS_Level3_Base, product_name='combine_dithers_all_exposures')\n", - "\n", - "asnallfile = 'for_spec3_all.json'\n", - "with open(asnallfile, 'w') as fpall:\n", - " fpall.write(asnall.dump()[1])" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "# Execute calwebb_spec3 pipeline. This will overwrite previous reductions.\n", - "\n", - "pipe3ss = Spec3Pipeline()\n", - "pipe3ss.master_background.skip = True\n", - "pipe3ss.mrs_imatch.skip = True\n", - "pipe3ss.outlier_detection.skip = True\n", - "pipe3ss.resample_spec.skip = True\n", - "pipe3ss.combine_1d.skip = True\n", - "pipe3ss.use_source_posn = 'True'\n", - "pipe3ss.subtract_background = 'True'\n", - "pipe3ss.output_file = 'allspec3'\n", - "pipe3ss.run(asnallfile)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Now to detect the point source in the detector image plane and extract and plot the spectra for each source\n", - "Steps involved:\n", - "1. Convert MJy/sr to Jy\n", - "2. Determine PSF centroid on the detector plane, if this one is not provided\n", - "3. Based on the across-slice position of the point source, determine the transmission and wavelength correction factors and apply correction\n", - "4. Project MRS PSF model from 3D spectral cube to 2D detector image plane\n", - "5. Determine pixel weights based on PSF and pixel variance information\n", - "6. Perform detector-based spectral extraction\n", - "7. If MRS spectral band is 3A (~12.2 micron), determine diffraction grating second order spectral leak contribution (from ~6.1 microns) and apply correction." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# centroid placeholders in local MRS coordinates\n", - "alpha_centers2D_dict = {}\n", - "beta_centers2D_dict = {}\n", - "# spectrum placeholder\n", - "isolambda_spec_optimal = {}" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# define spectral band\n", - "band = '1B'\n", - "ch = band[0]\n", - "\n", - "if ch in ['1', '2']:\n", - " det_id = 'MIRIFUSHORT'\n", - " band_combo = '12'\n", - "elif ch in ['3', '4']:\n", - " det_id = 'MIRIFULONG'\n", - " band_combo = '34'\n", - " \n", - "if band[1] == 'A':\n", - " band_id = 'short'\n", - " band_id_nr = 1\n", - "elif band[1] == 'B':\n", - " band_id = 'medium'\n", - " band_id_nr = 2\n", - "elif band[1] == 'C':\n", - " band_id = 'long'\n", - " band_id_nr = 3" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# load data file SCI and ERR extensions\n", - "if band in ['1A', '2A']:\n", - " hdu = fits.open('det_image_seq1_MIRIFUSHORT_12SHORTexp1120546_cal.fits')\n", - " fn = download_file('https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/jwst_miri_photom_0052.fits', cache=False)\n", - " hdu_photom = fits.open(fn) \n", - "elif band in ['1B', '2B']:\n", - " hdu = fits.open('det_image_seq1_MIRIFUSHORT_12MEDIUMexp1123047_cal.fits')\n", - " fn = download_file('https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/jwst_miri_photom_0064.fits', cache=False)\n", - " hdu_photom = fits.open(fn)\n", - "elif band in ['1C', '2C']:\n", - " hdu = fits.open('det_image_seq1_MIRIFUSHORT_12LONGexp1125354_cal.fits')\n", - " fn = download_file('https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/jwst_miri_photom_0060.fits', cache=False)\n", - " hdu_photom = fits.open(fn)\n", - "sci_img = hdu['SCI'].data\n", - "err_img = hdu['ERR'].data\n", - "\n", - "# load photom PIXSIZ extension to convert MJy/sr to Jansky\n", - "pixsiz_img = hdu_photom['PIXSIZ'].data\n", - "\n", - "# convert MJy/sr to Jansky for both the SCI and ERR extensions\n", - "sci_img = sci_img*pixsiz_img # MJy/sr --> MJy\n", - "err_img = err_img*pixsiz_img # MJy/sr --> MJy\n", - "\n", - "sci_img = sci_img*10**6 # MJy --> Jy\n", - "err_img = err_img*10**6 # MJy --> Jy\n", - "\n", - "# close opened fits files\n", - "hdu.close()\n", - "hdu_photom.close()\n", - "\n", - "# Point source corrections CDP values\n", - "# Transmission vs across-slice correction\n", - "fn = download_file('https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/miri-ifuptcor-59645.fits', cache=False)\n", - "hdu = fits.open(fn)\n", - "list_channels = ['CH1', 'CH2', 'CH3', 'CH4']\n", - "lamb_min = np.round(hdu['TRACOR'].data['WAVE_MIN'], 2)\n", - "lamb_max = np.round(hdu['TRACOR'].data['WAVE_MAX'], 2)\n", - "Tcen_min = np.round(hdu['TRACOR'].data['T_WMIN_CENTRE'], 2)\n", - "Tcen_max = np.round(hdu['TRACOR'].data['T_WMAX_CENTRE'], 2)\n", - "Tedg_min = np.round(hdu['TRACOR'].data['T_WMIN_EDGE'], 2)\n", - "Tedg_max = np.round(hdu['TRACOR'].data['T_WMAX_EDGE'], 2)\n", - "\n", - "# Wavelength offset vs across-slice correction\n", - "bands = hdu['WAVCORR_OPTICAL'].data['SUB_BAND']\n", - "x_slice_min_LT = np.round(hdu['WAVCORR_XSLICE'].data['XSLICE_MIN'][0], 8) # micron * arcsec^-1\n", - "x_slice_max_LT = np.round(hdu['WAVCORR_XSLICE'].data['XSLICE_MAX'][0], 8) # micron * arcsec^-1\n", - "\n", - "beta_slice = np.round(hdu['WAVCORR_OPTICAL'].data['BETA_SLICE'], 3)\n", - "\n", - "wave_min = np.round(hdu['WAVCORR_OPTICAL'].data['WAVE_MIN'], 2)\n", - "wave_max = np.round(hdu['WAVCORR_OPTICAL'].data['WAVE_MAX'], 2)\n", - "srp_min = hdu['WAVCORR_OPTICAL'].data['SRP_MIN']\n", - "srp_max = hdu['WAVCORR_OPTICAL'].data['SRP_MAX']\n", - "\n", - "beta_off = np.round(hdu['WAVCORR_SHIFT'].data['BETA_OFF'], 3)\n", - "Delta_s_min_LT = hdu['WAVCORR_SHIFT'].data['DS_MIN']\n", - "Delta_s_max_LT = hdu['WAVCORR_SHIFT'].data['DS_MAX']\n", - "\n", - "# Spectral leak (second order spectral response in form of an optical system transmission)\n", - "sys_transm_img = hdu['SCI'].data\n", - "hdu.close()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# define path to distortion (D2C) files (and create folder if it doesn't exist)\n", - "workDir = './miri_devel/'\n", - "cdpDir = workDir+'cdp_data/'\n", - "d2cDir = cdpDir+'DISTORTION/'\n", - "\n", - "# Evaluate largest distortion solution on the 2D detector (lowest slice transmission)\n", - "d2cMaps = d2cMapping(band, d2cDir, slice_transmission='10pc', fileversion=\"8B.05.02\")\n", - "sliceMap = d2cMaps['sliceMap']\n", - "lambdaMap = d2cMaps['lambdaMap']\n", - "alphaMap = d2cMaps['alphaMap']\n", - "betaMap = d2cMaps['betaMap']\n", - "nslices = d2cMaps['nslices']\n", - "\n", - "# some auxiliary data\n", - "det_dims = (1024, 1032)\n", - "bandlims = [lambdaMap[np.nonzero(lambdaMap)].min(), lambdaMap[np.nonzero(lambdaMap)].max()]\n", - "\n", - "# spectral grid based on which the point source centroid is determined\n", - "lambmin = mrs_aux(band)[3][0]\n", - "lambmax = mrs_aux(band)[3][1]\n", - "lambcens = np.arange(lambmin, lambmax, (lambmax-lambmin)/1024.)\n", - "lambfwhms = np.ones(len(lambcens))*(5*(lambmax-lambmin)/1024.)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Determine point source centroid\n", - "sign_amp2D, alpha_centers2D, beta_centers2D, sigma_alpha2D, sigma_beta2D, bkg_amp2D = point_source_centroiding(band, sci_img, d2cMaps, spec_grid=[lambcens, lambfwhms], fit='2D')\n", - "\n", - "alpha_centers2D_dict[band] = alpha_centers2D\n", - "beta_centers2D_dict[band] = beta_centers2D" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Diagnostic plot for the point source centroid: We care about the mean value of alpha, beta. These will be used for the PSF projection on the detector plane.\n", - "# Conceivably a PSF residual algorithm could be coded, to determine the alpha,beta centroid more accurately still.\n", - "fig, axs = plt.subplots(2, 1, figsize=(15, 10))\n", - "axs[0].plot(lambcens[20:-20], alpha_centers2D_dict[band][20:-20])\n", - "axs[0].plot(lambcens[20:-20], np.ones(len(lambcens[20:-20]))*np.mean(alpha_centers2D_dict[band][20:-20]), 'k--')\n", - "axs[1].plot(lambcens[20:-20], beta_centers2D_dict[band][20:-20])\n", - "axs[1].plot(lambcens[20:-20], np.ones(len(lambcens[20:-20]))*np.mean(beta_centers2D_dict[band][20:-20]), 'k--')\n", - "axs[0].set_ylabel(r'Along-slice position $\\alpha$ [arcsec]')\n", - "axs[0].set_title('MRS Band {}'.format(band))\n", - "axs[1].set_xlabel(r'Wavelength [$\\mu m$]')\n", - "axs[1].set_ylabel(r'Across-slice position $\\beta$ [arcsec]')\n", - "plt.tight_layout()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Determine transmission factor corresponding to source across-slice position $\\beta$, the signal should be multiplied by this factor" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "beta = np.mean(beta_centers2D_dict[band][20:-20])\n", - "wav_array = {}\n", - "Tb_min = {}\n", - "Tb_max = {}\n", - "Tdiff_point = {}\n", - "Tdiff_extended = {}\n", - "for i, key in enumerate(['CH1', 'CH2', 'CH3', 'CH4']):\n", - " wav_array[key] = np.arange(lamb_min[i], lamb_max[i], 0.01) # micron\n", - " Tb_min[key] = Tcen_min[i] + 2 * abs(beta) * (Tedg_min[i]-Tcen_min[i])\n", - " Tb_max[key] = Tcen_max[i] + 2 * abs(beta) * (Tedg_max[i]-Tcen_max[i])\n", - " \n", - " Tdiff_point[key] = Tb_min[key] + ((wav_array[key]-lamb_min[i])/(lamb_max[i]-lamb_min[i])) * (Tb_max[key]-Tb_min[key])\n", - " Tdiff_extended[key] = (Tcen_min[i]+Tedg_max[i])/2. + ((wav_array[key]-lamb_min[i])/(lamb_max[i]-lamb_min[i])) * ((Tcen_max[i]+Tedg_max[i]) - (Tcen_min[i]+Tedg_min[i])) / 2." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plt.figure(figsize=(12, 7))\n", - "for key in ['CH1', 'CH2', 'CH3', 'CH4']:\n", - " plt.plot(wav_array[key], Tdiff_point[key], 'k', linestyle='dashed')\n", - " plt.plot(wav_array[key], Tdiff_extended[key], 'k')\n", - "plt.xlabel(r'Wavelength [$\\mu m$]')\n", - "plt.ylabel(r'Transmission correction [%]')\n", - "legend_elements = [Line2D([0], [0], color='k', linestyle='dashed', label=r'Point source correction ($\\beta$ = {} arcsec)'.format(np.round(beta, 2))),\n", - " Line2D([0], [0], color='k', linestyle='solid', label='Extended source correction')]\n", - "plt.legend(handles=legend_elements, fontsize=20)\n", - "plt.tight_layout()\n", - "\n", - "plt.figure(figsize=(12, 7))\n", - "for key in ['CH1', 'CH2', 'CH3', 'CH4']:\n", - " plt.plot(wav_array[key], Tdiff_point[key]/Tdiff_extended[key], 'k', linestyle='dotted')\n", - "plt.xlabel(r'Wavelength [$\\mu m$]')\n", - "plt.ylabel(r'$\\zeta$ = Tdiff_point / Tdiff_extended')\n", - "legend_elements = [Line2D([0], [0], color='k', linestyle='dotted', label=r'$\\beta$ = {} arcsec'.format(np.round(beta, 2)))]\n", - "plt.legend(handles=legend_elements)\n", - "plt.tight_layout()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "key = 'CH' + band[0]\n", - "ip_zeta = interp1d(wav_array[key], Tdiff_point[key]/Tdiff_extended[key], kind='linear')\n", - "\n", - "zetaMap = np.ones(lambdaMap.shape)\n", - "sel = (sliceMap > 100*int(band[0])) & (sliceMap < 100*(int(band[0]) + 1))\n", - "zetaMap[sel] = ip_zeta(lambdaMap[sel])\n", - "\n", - "plt.figure(figsize=(8, 6))\n", - "plt.imshow(zetaMap, vmin=zetaMap[zetaMap != 0].min(), vmax=zetaMap[zetaMap != 0].max())\n", - "clb = plt.colorbar()\n", - "clb.set_label(r'$\\zeta$')\n", - "plt.xlabel('X-pixel')\n", - "plt.ylabel('Y-pixel')\n", - "plt.tight_layout()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Apply correction to science and error signal\n", - "sci_img *= zetaMap\n", - "err_img *= zetaMap" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Determine wavelength correction due to source across-slice position $\\beta$" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "beta_center = np.mean(beta_centers2D_dict[band][20:-20])\n", - "\n", - "# determine unique beta values for each slice\n", - "sel = np.where((sliceMap > 100*int(band[0])) & (sliceMap < 100*(int(band[0])+1)))\n", - "unique_betas = np.unique(betaMap[sel])\n", - "\n", - "# find index of value nearest to beta center in the unique_betas array\n", - "idx = np.abs(unique_betas-beta_center).argmin() \n", - "source_center_slice = idx+1\n", - "\n", - "# find index of value nearest to beta center in the unique_betas array\n", - "beta_brightest_slice = unique_betas[idx]\n", - "\n", - "# Point source offset from slice centre (-0.5 < tgtOffset < +0.5)\n", - "betaOffsetSlice0 = (beta_center-beta_brightest_slice)/d2cMaps['bdel'] \n", - "\n", - "# Define input data characteristics.\n", - "nSubBands = 3\n", - "optIndex = nSubBands * (int(band[0]) - 1) + (band_id_nr - 1)\n", - "betaSlice = beta_slice[optIndex]\n", - "\n", - "# Compute wavelength correction map (based on MIRI-DD-00006-ATC issue 4)\n", - "# The correction is valid only for the two slices nearest to the brightest slice\n", - "wavcorr_across_slice = np.zeros((1024, 1032))\n", - "for islice in range(-2, 3):\n", - " sel = (sliceMap == 100*int(band[0])+source_center_slice+islice)\n", - " tgtWave = lambdaMap[sel]\n", - " \n", - " # Step 1. Find target offset from slice centre in beta direction'\n", - " betaOff = betaOffsetSlice0 - islice\n", - "\n", - " # Step 2. Calculate Xslice\n", - " xSlice = betaSlice / tgtWave\n", - " xSliceref = betaSlice / ((lambmin+lambmax)/2.) # We define a reference xslice for the entire band to avoid jumps between beta_off values as a function of wavelength in Step 4.\n", - "\n", - " # Step 3. Derive scaled offset\n", - " betaOff_scaled = betaOff * xSlice\n", - " betaOff_scaled_ref = betaOff * xSliceref # We define a reference betaOff_scaled_ref for the entire band to avoid jumps between beta_off values as a function of wavelength in Step 4.\n", - "\n", - " # Step 4. Find look-up table betaOff_LT values which bracket betaOFF_scaled\n", - " nShiftValues = len(beta_off)\n", - " betaSign = 1\n", - " if (betaOff < 0.0):\n", - " betaSign = -1\n", - " betaOff_scaled = -1.0 * betaOff_scaled\n", - " indexA = -1 # Index in look up table_WCORR_SHIFT \n", - " for j in range(nShiftValues):\n", - " if (betaOff_scaled_ref < beta_off[j]):\n", - " if (indexA == -1):\n", - " indexA = j-1\n", - " indexB = indexA + 1\n", - " betaOff_LT_A = beta_off[indexA]\n", - " betaOff_LT_B = beta_off[indexB]\n", - "\n", - " # Step 5. Find ds_min and ds_max\n", - " betaFactor = (betaOff_scaled - betaOff_LT_A) / (betaOff_LT_B - betaOff_LT_A) \n", - "\n", - " ds_minA = Delta_s_min_LT[indexA]\n", - " ds_minB = Delta_s_min_LT[indexB]\n", - " ds_min = ds_minA + betaFactor * (ds_minB - ds_minA)\n", - " ds_min = betaSign * ds_min\n", - "\n", - " ds_maxA = Delta_s_max_LT[indexA]\n", - " ds_maxB = Delta_s_max_LT[indexA]\n", - " ds_max = ds_maxA + betaFactor * (ds_maxB - ds_maxA)\n", - " ds_max = betaSign * ds_max\n", - "\n", - " # Step 6. Find ds by interpolation between ds_min and ds_max\n", - " xSlice_min = x_slice_min_LT\n", - " xSlice_max = x_slice_max_LT\n", - " xFactor = (xSlice - xSlice_min) / (xSlice_max - xSlice_min)\n", - "\n", - " ds = ds_min + xFactor * (ds_max - ds_min)\n", - "\n", - " # Step 7. Convert ds to wavelength shift by intepolstion between SRP values\n", - " w_min = wave_min[optIndex]\n", - " w_max = wave_max[optIndex]\n", - " wFactor = (tgtWave - w_min) / (w_max - w_min)\n", - "\n", - " srp = srp_min[optIndex] + wFactor * (srp_max[optIndex] - srp_min[optIndex])\n", - "\n", - " wavcorr_across_slice[sel] = ds * tgtWave / srp\n", - "\n", - "# Plot wavelength correction map\n", - "plt.figure(figsize=(10, 6))\n", - "plt.imshow(wavcorr_across_slice)\n", - "clb = plt.colorbar()\n", - "clb.set_label('Wavelength correction [$\\u03bc m$]')\n", - "plt.xlabel('X-pixel')\n", - "plt.ylabel('Y-pixel')\n", - "plt.title('MRS spectral band {}'.format(band))\n", - "plt.tight_layout()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Apply correction to wavelength map\n", - "lambdaMap += wavcorr_across_slice" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Load MRS PSF model in the 3D spectral cube and project it on the 2D detector" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "psf_fits_file = cdpDir+\"MIRI_FM_{}_{}{}_PSF_07.02.00.fits\".format(det_id, ch, band_id.upper())\n", - "psffits = fits.open(psf_fits_file)\n", - "\n", - "# normalize and project PSF on 2D detector (normalization ensures that each cube slice of the model has a signal sum of 1)\n", - "psf_img = evaluate_psf_cdp(psffits, d2cMaps, source_center=[np.mean(alpha_centers2D_dict[band][20:-20]), np.mean(beta_centers2D_dict[band][20:-20])], norm=True)\n", - "\n", - "# close fits file\n", - "psffits.close()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Diagnostic plot to show how well the projected PSF matches the data\n", - "# For low SNR observations this will likely not be very useful\n", - "plt.close('all')\n", - "plt.figure(figsize=(16, 6))\n", - "if band[0] in ['1', '4']:\n", - " plt.plot(sci_img[512, 10:516], label='MIRISim signal')\n", - " scale_factor = (sci_img[512, 10:516][~np.isnan(sci_img[512, 10:516])].max()/psf_img[512, 10:516].max())\n", - " plt.plot(psf_img[512, 10:516]*scale_factor, label='PSF model (scaled)')\n", - "elif band[0] in ['2', '3']:\n", - " plt.plot(sci_img[512, 516:-10], label='MIRISim signal')\n", - " scale_factor = (sci_img[512, 516:-10][~np.isnan(sci_img[512, 516:-10])].max()/psf_img[512, 516:-10].max())\n", - " plt.plot(psf_img[512, 516:-10]*scale_factor, label='PSF model (scaled)')\n", - "plt.xlabel('X-pixel')\n", - "plt.ylabel('Signal [DN/sec]')\n", - "plt.legend()\n", - "plt.title('MRS spectral band {}'.format(band))\n", - "plt.tight_layout()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Determine optimal pixel weights for detector-based spectral extraction\n", - "weight_map = psf_img**2 / (err_img)**2" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Wavelength array used for the detector-based spectral extraction\n", - "# The grid step is approximately half the smallest spectral size in each band (due to the MRS distortion and slice curvature, different pixels contribute to each bin)\n", - "wav_array = {'1A': np.arange(4.9, 5.74, 0.0005), '1B': np.arange(5.67, 6.6, 0.0008), '1C': np.arange(6.45, 7.5, 0.001),\n", - " '2A': np.arange(7.477, 8.765, 0.0014), '2B': np.arange(8.711, 10.228, 0.0017), '2C': np.arange(10.017, 11.753, 0.002),\n", - " '3A': np.arange(11.481, 13.441, 0.0023), '3B': np.arange(13.319, 15.592, 0.0026), '3C': np.arange(15.4, 18.072, 0.0030),\n", - " '4A': np.arange(17.651, 20.938, 0.0036), '4B': np.arange(20.417, 24.22, 0.0042), '4C': np.arange(23.884, 28.329, 0.0048)}\n", - "nwavs = len(wav_array[band])\n", - "\n", - "# Evaluate smallest distortion solution on the 2D detector (largest slice transmission of 90%)\n", - "slice_transmission = '90pc'\n", - "d2cMaps = d2cMapping(band, d2cDir, slice_transmission=slice_transmission, fileversion=\"8B.05.02\")\n", - "lambdaMap = d2cMaps['lambdaMap']\n", - "lambdas = lambdaMap[np.nonzero(lambdaMap)].flatten()\n", - "\n", - "# Need approximate number of pixels contributing to each spectral bin\n", - "# This is because only one pixel is used for each detector column --> we want to avoid sampling issues\n", - "if band[0] in ['1', '2']:\n", - " if slice_transmission == '90pc':\n", - " npix = 380\n", - " elif slice_transmission == '10pc':\n", - " npix = 500\n", - "elif band[0] in ['3', '4']:\n", - " npix = 380\n", - "k = np.arange(npix)\n", - "\n", - "# Define pixel grid\n", - "X, Y = np.meshgrid(np.arange(1032), np.arange(1024))\n", - "X_flat = X[np.nonzero(lambdaMap)].flatten()\n", - "Y_flat = Y[np.nonzero(lambdaMap)].flatten()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Omit NaNs from analysis\n", - "nan_idx = ~np.isnan(sci_img)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Diffraction grating second order spectral leak contribution\n", - "if band == '1B':\n", - " \n", - " np.save('spectral_leak_transmission.npy', sys_transm_img)\n", - " \n", - " plt.close('all')\n", - " plt.figure(figsize=(8, 6))\n", - " plt.imshow(sys_transm_img*100.)\n", - " clb = plt.colorbar()\n", - " clb.set_label('Dichroic transmission [%]')\n", - " plt.xlabel('X-pixel')\n", - " plt.ylabel('Y-pixel')\n", - " plt.tight_layout()\n", - " \n", - " # Compute spectral leak signal from band 1B and update the error map\n", - " spectral_leak_sci_img = sci_img*sys_transm_img\n", - " spectral_leak_err_img = err_img*sys_transm_img\n", - " \n", - " plt.figure(figsize=(8, 6))\n", - " plt.imshow(spectral_leak_sci_img*1000)\n", - " clb = plt.colorbar()\n", - " clb.set_label('Spectral leak signal [mJy]')\n", - " plt.xlabel('X-pixel')\n", - " plt.ylabel('Y-pixel')\n", - " plt.tight_layout()\n", - " \n", - " # Redetermine optimal pixel weights for detector-based spectral extraction\n", - " spectral_leak_weight_map = psf_img**2 / (spectral_leak_err_img)**2\n", - " \n", - " # Initialize placeholder\n", - " spectral_leak_spectrum = np.zeros(nwavs)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Perform detector-based spectral extraction\n", - "isolambda_spec_optimal[band] = np.zeros(nwavs)\n", - "\n", - "for i in range(nwavs):\n", - " # determine pixels contributing to spectral bin\n", - " test_img = np.zeros((1024, 1032))\n", - " idxs = np.argsort(np.abs(lambdas-wav_array[band][i]))[k]\n", - " test_img[Y_flat[idxs], X_flat[idxs]] = 1\n", - " \n", - " # ascertain that only one pixel is selected per detector column\n", - " outliers = np.where(np.sum(test_img, axis=0) == 2)[0]\n", - "\n", - " for col in outliers:\n", - " Y_ = np.argsort(np.abs(lambdaMap[:, col] - wav_array[band][i]))[1]\n", - " test_img[Y_, col] = 0\n", - " \n", - " # plot meant for debugging !!!ONLY USE AT A SINGLE WAVELENGTH!!!\n", - " # plt.figure()\n", - " # plt.imshow(test_img,aspect=0.5)\n", - " # plt.tight_layout()\n", - "\n", - " # spectral extraction\n", - " sel_valid = (test_img == 1) & nan_idx\n", - " isolambda_spec_optimal[band][i] = np.sum(weight_map[sel_valid] * sci_img[sel_valid] / psf_img[sel_valid]) / np.sum(weight_map[sel_valid])\n", - " \n", - " if band == '1B':\n", - " # Estimate diffraction grating second order spectral leak contribution\n", - " spectral_leak_spectrum[i] = np.sum(spectral_leak_weight_map[sel_valid] * spectral_leak_sci_img[sel_valid] / psf_img[sel_valid]) / np.sum(spectral_leak_weight_map[sel_valid])" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# debugging spectral bins\n", - "\n", - "i = 200\n", - "# determine pixels contributing to spectral bin\n", - "test_img = np.zeros((1024, 1032))\n", - "idxs = np.argsort(np.abs(lambdas-wav_array[band][i]))[k]\n", - "test_img[Y_flat[idxs], X_flat[idxs]] = 1\n", - "\n", - "# ascertain that only one pixel is selected per detector column\n", - "outliers = np.where(np.sum(test_img, axis=0) == 2)[0]\n", - "\n", - "for col in outliers:\n", - " Y_ = np.argsort(np.abs(lambdaMap[:, col]-wav_array[band][i]))[1]\n", - " test_img[Y_, col] = 0\n", - "\n", - "plt.figure(figsize=(8, 6))\n", - "plt.imshow(test_img)\n", - "plt.imshow(d2cMaps['sliceMap'], alpha=0.4)\n", - "plt.xlabel('X-pixel')\n", - "plt.ylabel('Y-pixel')\n", - "plt.title('MRS spectral band {}'.format(band))\n", - "plt.tight_layout()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Compare resulting extracted spectrum to result from JWST-pipeline-reconstructed spectrum and original ISO spectrum\n", - "# Disclaimer: Due to flux conservation issues in MIRISim, the MIRISim MRS PSF is broader than the MRS PSF CDP model. \n", - "# This results on the detector-based spectrum having a significantly higher spectral baseline.\n", - "\n", - "# JWST pipeline result\n", - "hdu = fits.open('combine_dithers_all_exposures_ch{}-{}_x1d.fits'.format(ch, band_id))\n", - "wav = hdu[1].data['WAVELENGTH']\n", - "flux = hdu[1].data['FLUX']\n", - "hdu.close()\n", - "\n", - "# Original ISO spectrum\n", - "W_Per_hdu = fits.open('W_Per_spectrum_sws.fit')\n", - "wav_WPer = W_Per_hdu[0].data[:, 0]\n", - "flux_WPer = W_Per_hdu[0].data[:, 1]\n", - "sel_WPer = (wav_WPer > lambmin) & (wav_WPer < lambmax)\n", - "W_Per_hdu.close()\n", - "\n", - "# Plot\n", - "if band in ['1A', '1B', '1C']:\n", - " fudge_factor = 28.5\n", - "elif band in ['2A', '2B', '2C']:\n", - " fudge_factor = 45\n", - "plt.close('all')\n", - "plt.figure(figsize=(16, 6))\n", - "plt.plot(wav, flux, label='combine_dithers_all_exposures')\n", - "plt.plot(wav_array[band], isolambda_spec_optimal[band]/fudge_factor, label='Single dither detector-based extraction (scaled)')\n", - "plt.plot(wav_WPer[sel_WPer], flux_WPer[sel_WPer]/1000., label='ISO SWS spectrum of W Per')\n", - "plt.xlabel(r'Wavelength' u'[$\\u03bc m$]')\n", - "plt.ylabel(r'Signal [Jy]')\n", - "plt.title('MRS spectral band {}'.format(band))\n", - "plt.legend(fontsize=15)\n", - "plt.tight_layout()\n", - "\n", - "if band == '1B':\n", - " plt.figure(figsize=(16, 6))\n", - " plt.plot(wav_array[band], isolambda_spec_optimal[band], label='Original signal band 1B')\n", - " plt.plot(wav_array[band], spectral_leak_spectrum, label='Estimated spectral leak')\n", - " plt.xlabel('Wavelength [$\\u03bc m$]')\n", - " plt.ylabel('Spectral flux density [Jy]')\n", - " plt.title('MRS spectral band 1B')\n", - " plt.legend()\n", - " plt.tight_layout()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Correct for spectral leak if MRS band is set to 3A\n", - "## CAUTION: In order for the following cell to work, the user must have run this notebook already once for band 1B." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "if band == '3A':\n", - " spectral_leak_corr_band3A = isolambda_spec_optimal[band]-spectral_leak_spectrum" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# THE END\n", - "\n", - "## Still to do:\n", - "- Loop the notebook over all dither positions and average the resulting spectra. The reason why we may not want to use all four dithers in one go (i.e. to run the detector based extraction only once) is because we will likely find variations in the PSF as a function of which slice gets illuminated on the detector (think incoming cone angle and detector-scattered light).\n", - "- Consider updating current centroiding algorithm to use the cross-correlation method of the SDSS MANGA pipeline." - ] - }, - { - "cell_type": "markdown", - "metadata": { - "slideshow": { - "slide_type": "slide" - } - }, - "source": [ - "## About this notebook\n", - "**Author:** Ioannis (Yannis) Argyriou, Post-Doctoral Researcher, Institute of Astronomy, KU Leuven, Belgium \n", - "**Date:** 2022-01-08" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "***" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "[Top of Page](#top)" - ] - } - ], - "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.10" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/notebooks/MIRI/MRS_Mstar_analysis/JWST_Mstar_dataAnalysis_runpipeline.ipynb b/notebooks/MIRI/MRS_Mstar_analysis/JWST_Mstar_dataAnalysis_runpipeline.ipynb deleted file mode 100644 index f75c0e2b1..000000000 --- a/notebooks/MIRI/MRS_Mstar_analysis/JWST_Mstar_dataAnalysis_runpipeline.ipynb +++ /dev/null @@ -1,545 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "
\n", - "

This notebook is outdated and will be deprecated in October 2024

\n", - "
" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "slideshow": { - "slide_type": "slide" - } - }, - "source": [ - "# MRS Pipeline" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "**Use case:** Extract spatial-spectral features from IFU cube and measure their attributes.
\n", - "**Data:** Simulated [MIRI MRS](https://jwst-docs.stsci.edu/mid-infrared-instrument/miri-observing-modes/miri-medium-resolution-spectroscopy) spectrum of AGB star.
\n", - "**Tools:** specutils, jwst, photutils, astropy, aplpy, scipy.
\n", - "**Cross-intrument:** NIRSpec, MIRI.
\n", - "**Documentation:** This notebook is part of a STScI's larger [post-pipeline Data Analysis Tools Ecosystem](https://jwst-docs.stsci.edu/jwst-post-pipeline-data-analysis) and can be [downloaded](https://github.com/spacetelescope/dat_pyinthesky/tree/main/jdat_notebooks/MRS_Mstar_analysis) directly from the [JDAT Notebook Github directory](https://github.com/spacetelescope/jdat_notebooks).
\n", - "**Source of Simulations:** [MIRISim](https://www.stsci.edu/jwst/science-planning/proposal-planning-toolbox/mirisim)
\n", - "**Pipeline Version:** [JWST Pipeline](https://jwst-docs.stsci.edu/jwst-data-reduction-pipeline)
\n", - "\n", - "**Note**: This notebook includes MIRI simulated data cubes obtained using MIRISim (https://wiki.miricle.org//bin/view/Public/MIRISim_Public)\n", - "and run through the JWST pipeline (https://jwst-pipeline.readthedocs.io/en/latest/) of\n", - "point sources with spectra representative of late M type stars.\n", - "\n", - "## Introduction\n", - "\n", - "This notebook analyzes one star represented by a dusty SED corresponding to the ISO SWS spectrum of\n", - "W Per from Kraemer et al. (2002) and Sloan et al. (2003) to cover the MRS spectral range 5-28 microns. Analysis of JWST spectral cubes requires extracting spatial-spectral features of interest and measuring their attributes.\n", - "\n", - "This is the first notebook, which will process the data and automatically detect and extract spectra for the point source. The workflow will use `photutils` to automatically detect sources in the cube to extract the spectrum of the point source. Then it will read in the spectra generated at Stage 3 of the JWST pipeline.\n" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "slideshow": { - "slide_type": "slide" - } - }, - "source": [ - "## Import Packages" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "slideshow": { - "slide_type": "fragment" - } - }, - "outputs": [], - "source": [ - "# Import useful python packages\n", - "import numpy as np\n", - "\n", - "# Import packages to display images inline in the notebook\n", - "import matplotlib.pyplot as plt \n", - "%matplotlib inline \n", - "\n", - "# Set general plotting options\n", - "params = {'legend.fontsize': '18', 'axes.labelsize': '18', \n", - " 'axes.titlesize': '18', 'xtick.labelsize': '18', \n", - " 'ytick.labelsize': '18', 'lines.linewidth': 2, 'axes.linewidth': 2, 'animation.html': 'html5'}\n", - "plt.rcParams.update(params)\n", - "plt.rcParams.update({'figure.max_open_warning': 0})" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import warnings\n", - "warnings.simplefilter('ignore')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false, - "jupyter": { - "outputs_hidden": false - }, - "pycharm": { - "name": "#%%\n" - } - }, - "outputs": [], - "source": [ - "# Import astropy packages \n", - "from astropy import units as u\n", - "from astropy.io import ascii\n", - "from astropy.wcs import WCS\n", - "from astropy.table import Table, vstack\n", - "from astropy.stats import sigma_clipped_stats\n", - "from astropy.nddata import StdDevUncertainty\n", - "from astropy.io import fits # added by BAS on 8 April 2021\n", - "from astropy.utils.data import get_pkg_data_filename\n", - "\n", - "# To find stars in the MRS spectralcubes and do aperture photometry\n", - "from photutils.aperture import CircularAperture\n", - "from photutils.detection import DAOStarFinder\n", - "\n", - "# To deal with 1D spectrum\n", - "from specutils import Spectrum1D\n", - "from specutils.fitting import fit_generic_continuum\n", - "from specutils.manipulation import box_smooth, extract_region, SplineInterpolatedResampler\n", - "from specutils.analysis import line_flux, centroid, equivalent_width\n", - "from specutils.spectra import SpectralRegion\n", - "from specutils import SpectrumList\n", - "\n", - "# To fit a curve to the data\n", - "from scipy.optimize import curve_fit\n", - "\n", - "import os" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Set paths to the Data and Outputs\n", - "\n", - "Use MIRISim JWST pipeline processed data in future iterations." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "os.environ['CRDS_PATH'] = os.environ['HOME']+'/crds_cache'\n", - "os.environ['CRDS_SERVER_URL'] = 'https://jwst-crds-pub.stsci.edu'" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# import pipeline\n", - "\n", - "from jwst.pipeline import Detector1Pipeline\n", - "from jwst.pipeline import Spec2Pipeline\n", - "from jwst.pipeline import Spec3Pipeline\n", - "from jwst.extract_1d import Extract1dStep\n", - "import json\n", - "import glob\n", - "from jwst.associations.lib.rules_level3_base import DMS_Level3_Base\n", - "from jwst.associations import asn_from_list\n", - "import crds\n", - "from jdaviz.app import Application\n", - "import asdf\n", - "from photutils import aperture_photometry\n", - "import os" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Download data if you don't already have it.\n", - "\n", - "import urllib.request\n", - "\n", - "if os.path.exists(\"20210413_120546_mirisim.tar.gz\"):\n", - " print(\"20210413_120546_mirisim.tar.gz Exists\")\n", - "else:\n", - " url = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/20210413_120546_mirisim.tar.gz'\n", - " urllib.request.urlretrieve(url, './20210413_120546_mirisim.tar.gz')\n", - " \n", - "if os.path.exists(\"20210413_123047_mirisim.tar.gz\"):\n", - " print(\"20210413_123047_mirisim.tar.gz Exists\")\n", - "else:\n", - " url = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/20210413_123047_mirisim.tar.gz'\n", - " urllib.request.urlretrieve(url, './20210413_123047_mirisim.tar.gz')\n", - " \n", - "if os.path.exists(\"20210413_125354_mirisim.tar.gz\"):\n", - " print(\"20210413_125354_mirisim.tar.gz Exists\")\n", - "else:\n", - " url = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/20210413_125354_mirisim.tar.gz'\n", - " urllib.request.urlretrieve(url, './20210413_125354_mirisim.tar.gz')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Unzip Tar Files\n", - "\n", - "import tarfile\n", - "\n", - "# Unzip files if they haven't already been unzipped\n", - "if os.path.exists(\"20210413_120546_mirisim/\"):\n", - " print(\"20210413_120546_mirisim Exists\")\n", - "else:\n", - " tar = tarfile.open('./20210413_120546_mirisim.tar.gz', \"r:gz\")\n", - " tar.extractall()\n", - " tar.close()\n", - " \n", - "if os.path.exists(\"20210413_123047_mirisim/\"):\n", - " print(\"20210413_123047_mirisim Exists\")\n", - "else:\n", - " tar = tarfile.open('./20210413_123047_mirisim.tar.gz', \"r:gz\")\n", - " tar.extractall()\n", - " tar.close()\n", - " \n", - "if os.path.exists(\"20210413_125354_mirisim/\"):\n", - " print(\"20210413_125354_mirisim Exists\")\n", - "else:\n", - " tar = tarfile.open('./20210413_125354_mirisim.tar.gz', \"r:gz\")\n", - " tar.extractall()\n", - " tar.close()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Run Pipeline\n", - "\n", - "The various [stages of the pipeline](https://jwst-pipeline.readthedocs.io/en/latest/jwst/pipeline/main.html#pipelines) can be [run within Python](https://jwst-pipeline.readthedocs.io/en/latest/jwst/introduction.html#running-from-within-python). For a more in depth tutorial on running the pipelines, check out the [JWebbinars](https://www.stsci.edu/jwst/science-execution/jwebbinars)." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Execute calwebb_detector1 pipeline on raw simulation output. This will overwrite previous reductions.\n", - "\n", - "allshortfiles = glob.glob('20210413_*_mirisim/det_images/*MIRIFUSHORT*fits')\n", - "alllongfiles = glob.glob('20210413_*_mirisim/det_images/*MIRIFULONG*fits')\n", - " \n", - "pipe1short = Detector1Pipeline()\n", - "\n", - "# run calwebb_detector1 on the MIRIFUSHORT data separate from MIRIFULONG data, as it saves time this way\n", - "for shortfile in allshortfiles:\n", - " print(shortfile)\n", - " baseshort, remaindershort = shortfile.split('.')\n", - " \n", - " # If you run your own simulations, you will need to update these hardcoded files.\n", - " beforestuffshort, dateafterstuffshort = shortfile.split('20210413_') \n", - " datestringshort, afterstuffshort = dateafterstuffshort.split('_mirisim')\n", - " \n", - " pipe1short.refpix.skip = True\n", - " pipe1short.output_file = baseshort + datestringshort\n", - " \n", - " pipe1short.run(shortfile)\n", - "\n", - "pipe1long = Detector1Pipeline()\n", - "\n", - "for longfile in alllongfiles:\n", - " print(longfile)\n", - " baselong, remainderlong = longfile.split('.')\n", - " \n", - " # If you run your own simulations, you will need to update these hardcoded files.\n", - " beforestufflong, dateafterstufflong = longfile.split('20210413_')\n", - " datestringlong, afterstufflong = dateafterstufflong.split('_mirisim')\n", - " \n", - " pipe1long.refpix.skip = True\n", - " pipe1long.output_file = baselong + datestringlong\n", - " \n", - " pipe1long.run(longfile)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Execute calwebb_spec2 pipeline. This will overwrite previous reductions.\n", - "\n", - "# All the local calwebb_detector1 files\n", - "allshortfiles2 = glob.glob('det_image_*_MIRIFUSHORT_*_rate.fits')\n", - "alllongfiles2 = glob.glob('det_image_*_MIRIFULONG_*_rate.fits')\n", - "\n", - "for short2file in allshortfiles2:\n", - " print(short2file)\n", - " pipe2short = Spec2Pipeline()\n", - " base2short, remainder2short = short2file.split('.')\n", - " \n", - " pipe2short.straylight.skip = True\n", - " \n", - " # If you run your own simulations, you will need to update these hardcoded files.\n", - " if (short2file == 'det_image_seq1_MIRIFUSHORT_12LONGexp1125354_rate.fits'):\n", - " print('this one will have the level 2b cube built')\n", - " else:\n", - " pipe2short.cube_build.skip = True\n", - " pipe2short.extract_1d.skip = True\n", - " pipe2short.output_file = base2short\n", - " \n", - " pipe2short.run(short2file)\n", - "\n", - "for long2file in alllongfiles2:\n", - " print(long2file)\n", - " pipe2long = Spec2Pipeline()\n", - " base2long, remainder2long = long2file.split('.')\n", - " \n", - " pipe2long.straylight.skip = True\n", - " # If you run your own simulations, you will need to update these hardcoded files.\n", - " if (long2file == 'det_image_seq1_MIRIFULONG_34SHORTexp1120546_rate.fits'):\n", - " print('this one will have the level 2b cube built')\n", - " else:\n", - " pipe2long.cube_build.skip = True\n", - " pipe2long.extract_1d.skip = True\n", - " pipe2long.output_file = base2long\n", - " \n", - " pipe2long.run(long2file)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Now to detect the point source in the datacube and extract and plot the spectra for each source\n", - "\n", - "For data cubes like the JWST/MIRI MRS information on the point sources in the FOV and also obtaining a source subtracted\n", - " data cube will be necessary (See the `PampelMuse` software for an example on how spectral extraction is implemented for\n", - " near-IR data cubes like MUSE).\n", - "\n", - "Note these backgrounds of diffuse emission can be quite complex.\n", - "\n", - "On these source extracted data cubes (see `SUBTRES` in `PampelMuse`) I would like to produce moment maps\n", - "(https://casa.nrao.edu/Release3.4.0/docs/UserMan/UserManse41.html) and Position-Velocity (PV) diagrams\n", - "(https://casa.nrao.edu/Release4.1.0/doc/UserMan/UserManse42.html)." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### 1) Use `Photutils` to detect stars/point sources in the continuum image\n", - "\n", - "The first step of the analysis is to identify those sources for which it is feasible to extract spectra from the IFU\n", - "data. Ideally we can estimate the signal-to-noise ratio (S/N) for all sources in the cube, do a number of checks to\n", - "determine the status of every source and loop through these (brightest first) to extract the spectra. Open up the Level 2 Cubes and use photutils to search for point sources for Level 3 extraction." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# If you run your own simulations, you will need to update these hardcoded files.\n", - "l_cube_file = 'det_image_seq1_MIRIFULONG_34SHORTexp1120546_s3d.fits'\n", - "s_cube_file = 'det_image_seq1_MIRIFUSHORT_12LONGexp1125354_s3d.fits'\n", - "\n", - "with fits.open(s_cube_file) as hdu_s_cube:\n", - " s_cube = hdu_s_cube['SCI'].data\n", - " s_med_cube = np.nanmedian(s_cube, axis=0)\n", - "\n", - "mean, median, std = sigma_clipped_stats(s_med_cube, sigma = 2.0)\n", - "\n", - "# Get a list of sources using a dedicated source detection algorithm\n", - "# Find sources at least 3* background (typically)\n", - "\n", - "daofind = DAOStarFinder(fwhm = 2.0, threshold = 3. * std)\n", - "sources = daofind(s_med_cube - median) \n", - "print(\"\\n Number of sources in field: \", len(sources))\n", - "\n", - "# Positions in pixels\n", - "positions = Table([sources['xcentroid'], sources['ycentroid']])\n", - "\n", - "# Convert to RA & Dec (ICRS)\n", - "peakpixval = np.zeros(len(sources['xcentroid']))\n", - "for count_s, _ in enumerate(sources):\n", - " peakpixval[count_s] = s_med_cube[int(np.round(sources['xcentroid'][count_s])), int(np.round(sources['ycentroid'][count_s]))]\n", - "print('peak pixel x =')\n", - "print(sources['xcentroid'][np.argmax(peakpixval)])\n", - "print('peak pixel y =')\n", - "print(sources['ycentroid'][np.argmax(peakpixval)])\n", - "\n", - "plt.imshow(s_med_cube, vmin=0, vmax=100)#.value)\n", - "plt.tight_layout()\n", - "plt.scatter(sources['xcentroid'], sources['ycentroid'], c = \"red\", marker = \"+\", s=50)\n", - "plt.scatter(sources['xcentroid'][np.argmax(peakpixval)], sources['ycentroid'][np.argmax(peakpixval)], c = 'black', marker='+', s=50)\n", - "plt.show()\n", - "\n", - "f0 = fits.open(s_cube_file)\n", - "w0 = WCS(f0[('sci',1)].header, f0)\n", - "f0.close()\n", - "\n", - "radec = w0.all_pix2world([sources['xcentroid'][np.argmax(peakpixval)]], [sources['ycentroid'][np.argmax(peakpixval)]], [1], 1)\n", - "\n", - "# Take the brightest source flux and take that to be your primary point source for extraction\n", - "ra_ptsrc = radec[0][0]\n", - "dec_ptsrc = radec[1][0]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Due to the way the pipeline currently extracts Level3 data, you must update the headers to be centered on the point source of your choosing from the step above.\n", - "all_files = glob.glob('det_image_*_cal.fits')\n", - "targra = ra_ptsrc\n", - "targdec = dec_ptsrc\n", - "for thisfile in all_files:\n", - " base, remainder = thisfile.split('.')\n", - " outfilename = base + '_fix.' + remainder\n", - " print(outfilename)\n", - " \n", - " with fits.open(thisfile) as hduthis:\n", - " hduthis['SCI'].header['SRCTYPE'] = 'POINT'\n", - " hduthis[0].header['TARG_RA'] = targra\n", - " hduthis[0].header['TARG_DEC'] = targdec\n", - " hduthis.writeto(outfilename, overwrite=True)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# set up needed reference file(s) for spec3\n", - "\n", - "file_all_list = glob.glob('det_image_*_cal_fix.fits')\n", - "\n", - "asnall = asn_from_list.asn_from_list(file_all_list, rule=DMS_Level3_Base, product_name='combine_dithers_all_exposures')\n", - "\n", - "asnallfile = 'for_spec3_all.json'\n", - "with open(asnallfile, 'w') as fpall:\n", - " fpall.write(asnall.dump()[1])" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Execute calwebb_spec3 pipeline. This will overwrite previous reductions.\n", - "\n", - "pipe3ss = Spec3Pipeline()\n", - "pipe3ss.master_background.skip = True\n", - "pipe3ss.mrs_imatch.skip = True\n", - "pipe3ss.outlier_detection.skip = True\n", - "pipe3ss.resample_spec.skip = True\n", - "pipe3ss.combine_1d.skip = True\n", - "pipe3ss.use_source_posn = 'True'\n", - "pipe3ss.subtract_background = 'True'\n", - "pipe3ss.output_file = 'allspec3'\n", - "pipe3ss.run(asnallfile)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Next Step\n", - "\n", - "Proceed to Notebook 2 for visualization and data anlysis." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Additional Resources\n", - "\n", - "- [PampelMuse](https://gitlab.gwdg.de/skamann/pampelmuse)\n", - "- [CASA](https://casa.nrao.edu/Release3.4.0/docs/UserMan/UserManse41.html)" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "slideshow": { - "slide_type": "slide" - } - }, - "source": [ - "## About this notebook\n", - "**Author:** Olivia Jones, Project Scientist, UK ATC.\n", - "**Updated On:** 2020-08-11\n", - "**Later Updated On:** 2021-09-06 by B. Sargent, STScI Scientist, Space Telescope Science Institute" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "***" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "[Top of Page](#top)" - ] - } - ], - "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.10" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -}