Skip to content

Commit

Permalink
path fixes, print statements fixes, and adding helper functions
Browse files Browse the repository at this point in the history
  • Loading branch information
@Jackie-Brown
Jackie-Brown committed Oct 19, 2024
1 parent 321c704 commit 30c6de1
Showing 1 changed file with 66 additions and 55 deletions.
121 changes: 66 additions & 55 deletions notebooks/STIS/low_count_uncertainties/Low_Count_Uncertainties.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -96,6 +96,9 @@
"from astropy.stats import poisson_conf_interval\n",
"from scipy.stats import norm, poisson\n",
"\n",
"# Import for shortened outputs\n",
"from IPython.display import clear_output\n",
"\n",
"# Using a colorblind friendly style\n",
"plt.style.use('tableau-colorblind10')"
]
Expand Down Expand Up @@ -127,9 +130,14 @@
"# get a list of files assiciated with that target\n",
"target_list = Observations.get_product_list(target)\n",
"\n",
"# Download and read in only the x1d file\n",
"# download only the x1d file\n",
"Observations.download_products(target_list, extension='x1d.fits')\n",
"with fits.open('mastDownload/HST/obgh07020/obgh07020_x1d.fits') as hdu:\n",
"\n",
"# get file path\n",
"file_path_obgh07020 = os.path.join('mastDownload', 'HST', target_id, f'{target_id}_x1d.fits')\n",
"\n",
"# read in x1d\n",
"with fits.open(file_path_obgh07020) as hdu:\n",
" header = hdu[0].header\n",
" data = hdu[1].data\n",
" exptime = hdu[1].header['EXPTIME']\n",
Expand All @@ -147,11 +155,10 @@
"print(hdu[0].header['TELESCOP'], hdu[0].header['INSTRUME'],\n",
" hdu[0].header['DETECTOR'], hdu[0].header['OPT_ELEM'],\n",
" hdu[0].header['APERTURE'])\n",
"print('Program: ', hdu[0].header['PROPOSID'])\n",
"print('PI: ', hdu[0].header['PR_INV_F'], hdu[0].header['PR_INV_M'],\n",
" hdu[0].header['PR_INV_L'])\n",
"print('Target: ', hdu[0].header['TARGNAME'])\n",
"print('Time: ', hdu[0].header['TDATEOBS'], hdu[0].header['TTIMEOBS'])"
"print(f\"Program: {hdu[0].header['PROPOSID']}\")\n",
"print(f\"PI: {hdu[0].header['PR_INV_F']} {hdu[0].header['PR_INV_M']} {hdu[0].header['PR_INV_L']}\")\n",
"print(f\"Target: {hdu[0].header['TARGNAME']}\")\n",
"print(f\"Time: {hdu[0].header['TDATEOBS']} {hdu[0].header['TTIMEOBS']}\")"
]
},
{
Expand Down Expand Up @@ -534,9 +541,14 @@
"# get a list of files assiciated with that target\n",
"target_list = Observations.get_product_list(target)\n",
"\n",
"# Download and read in just the x1d file\n",
"# Download just the x1d file\n",
"Observations.download_products(target_list, extension='ld9m12erq_x1d.fits')\n",
"with fits.open('mastDownload/HST/ld9m12erq/ld9m12erq_x1d.fits') as hdu:\n",
"\n",
"# get file path\n",
"file_path_ld9m12erq = os.path.join('mastDownload', 'HST', 'ld9m12erq', 'ld9m12erq_x1d.fits')\n",
"\n",
"# read in x1d\n",
"with fits.open(file_path_ld9m12erq) as hdu:\n",
" header = hdu[0].header\n",
" data = hdu[1].data\n",
" exptime = hdu[1].header['EXPTIME']\n",
Expand All @@ -558,11 +570,10 @@
"print(hdu[0].header['TELESCOP'], hdu[0].header['INSTRUME'],\n",
" hdu[0].header['DETECTOR'], hdu[0].header['OPT_ELEM'],\n",
" hdu[0].header['APERTURE'])\n",
"print('Program: ', hdu[0].header['PROPOSID'])\n",
"print('PI: ', hdu[0].header['PR_INV_F'], hdu[0].header['PR_INV_M'],\n",
" hdu[0].header['PR_INV_L'])\n",
"print('Target: ', hdu[0].header['TARGNAME'])\n",
"print('Date: ', hdu[0].header['DATE'])"
"print(f\"Program: {hdu[0].header['PROPOSID']}\")\n",
"print(f\"PI: {hdu[0].header['PR_INV_F']} {hdu[0].header['PR_INV_M']} {hdu[0].header['PR_INV_L']}\")\n",
"print(f\"Target: {hdu[0].header['TARGNAME']}\")\n",
"print(f\"Date: {hdu[0].header['DATE']}\")"
]
},
{
Expand Down Expand Up @@ -738,10 +749,12 @@
"# Download the tag file for this observation\n",
"Observations.download_products(target_list, extension='tag.fits')\n",
"\n",
"# get file paths\n",
"tag_file = os.path.join('mastDownload', 'HST', target_id, f'{target_id}_tag.fits')\n",
"raw_file = os.path.join('mastDownload', 'HST', target_id, f'{target_id}_tag.fits')\n",
"\n",
"# Split up the time-tag file into 10 subexposures\n",
"stistools.inttag.inttag(\"./mastDownload/HST/obgh07020/obgh07020_tag.fits\",\n",
" \"./mastDownload/HST/obgh07020/obgh07020_raw.fits\",\n",
" rcount=10)\n",
"stistools.inttag.inttag(tag_file, raw_file, rcount=10)\n",
"\n",
"# We then want to re-reduced the raw data to get x1d files,\n",
"# and then we'll compare:\n",
Expand All @@ -753,19 +766,21 @@
"os.environ[\"CRDS_PATH\"] = crds_path\n",
"os.environ[\"CRDS_SERVER_URL\"] = \"https://hst-crds.stsci.edu\"\n",
"os.environ[\"oref\"] = os.path.join(crds_path, \"references/hst/oref/\")\n",
"!crds bestrefs --update-bestrefs --sync-references=1 --files ./mastDownload/HST/obgh07020/obgh07020_raw.fits\n",
"!crds bestrefs --update-bestrefs --sync-references=1 --files raw_file\n",
"\n",
"# We need to remove the x1d file we downloaded before\n",
"# because we'll replace it with the inttag version\n",
"inttag_path = os.path.join('mastDownload', 'HST', target_id, 'inttag*.fits')\n",
"try:\n",
" os.remove('./mastDownload/HST/obgh07020/inttag*.fits')\n",
" os.remove(inttag_path)\n",
"except OSError:\n",
" print('No need to remove...')\n",
"# Now run calstis...\n",
"wave_file = './mastDownload/HST/obgh07020/obgh07020_wav.fits'\n",
"raw_file = './mastDownload/HST/obgh07020/obgh07020_raw.fits'\n",
"wave_file = os.path.join('mastDownload', 'HST', target_id, f'{target_id}_wav.fits')\n",
"stistools.calstis.calstis(raw_file, wavecal=wave_file, verbose=False,\n",
" outroot=\"./mastDownload/HST/obgh07020/inttag\")"
" outroot=os.path.join('mastDownload', 'HST', target_id, 'inttag'))\n",
"\n",
"clear_output()"
]
},
{
Expand All @@ -783,8 +798,8 @@
"metadata": {},
"outputs": [],
"source": [
"hdu = fits.open('./mastDownload/HST/obgh07020/obgh07020_x1d.fits')\n",
"hdu_int = fits.open('./mastDownload/HST/obgh07020/inttag_x1d.fits')\n",
"hdu = fits.open(file_path_obgh07020)\n",
"hdu_int = fits.open(os.path.join('mastDownload', 'HST', target_id, 'inttag_x1d.fits'))\n",
"\n",
"plt.figure()\n",
"plt.plot(hdu[1].data['WAVELENGTH'][0],\n",
Expand Down Expand Up @@ -834,47 +849,43 @@
"scenario1 = [0, 0, 10, 0, 0]\n",
"scenario2 = [2, 2, 2, 2, 2]\n",
"\n",
"# define helper functions for the following print statements\n",
"def quadrature(scenario):\n",
" return np.sqrt(np.sum([np.sqrt(i)**2 for i in scenario]))\n",
"\n",
"def sum_sq(scenario):\n",
" return np.sqrt(np.sum(scenario))\n",
"\n",
"def pixel_err(scenario):\n",
" return np.sqrt(np.sum(\n",
" [(poisson_conf_interval(i, interval='frequentist-confidence')[1]-i)**2 \n",
" for i in scenario]))\n",
"\n",
"def upper_err(scenario):\n",
" return (poisson_conf_interval(\n",
" np.sum(scenario), interval='frequentist-confidence')[1]\n",
" - np.sum(scenario))\n",
"\n",
"print('Scenario 1: ', scenario1)\n",
"print('Scenario 2: ', scenario2)\n",
"\n",
"print('---')\n",
"print('Root-N Approximation')\n",
"print('---')\n",
"print('Scenario 1')\n",
"print('---')\n",
"print('Root-N approximation, add in quadrature: ',\n",
" np.sqrt(np.sum([np.sqrt(i)**2 for i in scenario1])))\n",
"print('Root-N approximation, sum column, then take sqrt: ',\n",
" np.sqrt(np.sum(scenario1)))\n",
"print('Root-N approximation, add in quadrature: ', quadrature(scenario1))\n",
"print('Root-N approximation, sum column, then take sqrt: ', sum_sq(scenario1))\n",
"print('---')\n",
"tmp = np.sqrt(\n",
" np.sum(\n",
" [(poisson_conf_interval(i, interval='frequentist-confidence')[1]-i)**2\n",
" for i in scenario1]))\n",
"print('Poisson, pixel-by-pixel error calculation: ', tmp)\n",
"\n",
"print('Poisson upper limit, sum column, then calculate error: ',\n",
" poisson_conf_interval(np.sum(scenario1),\n",
" interval='frequentist-confidence')[1] -\n",
" np.sum(scenario1))\n",
"print('Poisson, pixel-by-pixel error calculation: ', pixel_err(scenario1))\n",
"print('Poisson upper limit, sum column, then calculate error: ', upper_err(scenario1))\n",
"print('---')\n",
"print('Scenario 2')\n",
"print('---')\n",
"print('Root-N approximation, add in quadrature: ',\n",
" np.sqrt(np.sum([np.sqrt(i)**2 for i in scenario2])))\n",
"print('Root-N approximation, sum column, then take sqrt: ',\n",
" np.sqrt(np.sum(scenario2)))\n",
"print('Root-N approximation, add in quadrature: ', quadrature(scenario2))\n",
"print('Root-N approximation, sum column, then take sqrt: ', sum_sq(scenario2))\n",
"print('---')\n",
"print('Poisson, pixel-by-pixel error calculation: ',\n",
" np.sqrt(\n",
" np.sum(\n",
" [(poisson_conf_interval(\n",
" i, interval='frequentist-confidence')[1]-i)**2\n",
" for i in scenario2])))\n",
"print('Poisson upper limit, sum column, then calculate error: ',\n",
" poisson_conf_interval(\n",
" np.sum(scenario2), interval='frequentist-confidence')[1]\n",
" - np.sum(scenario2))"
"print('Poisson, pixel-by-pixel error calculation: ', pixel_err(scenario2))\n",
"print('Poisson upper limit, sum column, then calculate error: ', upper_err(scenario2))"
]
},
{
Expand Down Expand Up @@ -920,7 +931,7 @@
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"display_name": "stis",
"language": "python",
"name": "python3"
},
Expand All @@ -934,7 +945,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.6"
"version": "3.7.10"
}
},
"nbformat": 4,
Expand Down

0 comments on commit 30c6de1

Please sign in to comment.