diff --git a/notebooks/effective_survey_speed_trending.ipynb b/notebooks/effective_survey_speed_trending.ipynb new file mode 100644 index 0000000..0dd2740 --- /dev/null +++ b/notebooks/effective_survey_speed_trending.ipynb @@ -0,0 +1,675 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0fe4b48e-ede6-42d2-a1bd-21fe3eee0c15", + "metadata": {}, + "source": [ + "* slew performance\n", + "* delay to on-sky\n", + "* loss of time at end of night\n", + "* faults within the night\n", + "* weather downtime\n", + "* scheduled downtime\n", + "\n", + "Relating these to overall system availability:\n", + "\n", + "For dayobs 20251108\n", + "* total time in the night = 8.06 hours\n", + "* time lost from -12 deg twilight to start of science = 1.2 hours (includes BLOCK-T539)\n", + "* time lost due to faults and recovery ~2.44 hours\n", + "* time lost current slew speed compared to baseline_v5.1.0 = 0.5 hours\n", + "* time slewing/changing filter = 1.4 hours\n", + "* time that would be spent slewing/changing filter if slew like v5.1. = 0.87 hours\n", + "* time spent taking image = 2.9 hours\n", + "* time exposing at sky = 2.82 (0.094 hours shutter movement)\n", + "* science open shutter fraction 20251108 = 35%\n", + "* median open shutter fraction in v5.1 = 75%" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2d4b4c71-a0d1-4dfd-a243-e74568428213", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import pandas as pd\n", + "\n", + "import matplotlib.pyplot as plt\n", + "\n", + "%matplotlib widget" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "48a9111f-7037-4f88-ac09-aee75d506b15", + "metadata": {}, + "outputs": [], + "source": [ + "from rubin_scheduler.site_models import Almanac\n", + "from rubin_scheduler.utils import Site" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a8cf1335-dfcc-41e4-beab-c89c188973de", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "os.environ[\"no_proxy\"] += \",.consdb\"\n", + "\n", + "from lsst.summit.utils import ConsDbClient\n", + "\n", + "client = ConsDbClient(\"http://consdb-pq.consdb:8080/consdb\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8a3aff13-ae57-4613-87f2-1f062e1e10c7", + "metadata": {}, + "outputs": [], + "source": [ + "instrument = 'lsstcam'" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ded6d1f7-6801-49f1-8691-5f2397a29e7a", + "metadata": {}, + "outputs": [], + "source": [ + "def fetch(day_obs_min, day_obs_max):\n", + " visits_query = f'''\n", + " SELECT \n", + " *\n", + " FROM\n", + " cdb_{instrument}.visit1 AS visit1 \n", + " WHERE\n", + " airmass > 0.\n", + " AND img_type in ('acq', 'science', 'cwfs', 'focus')\n", + " AND day_obs >= {day_obs_min} AND day_obs <= {day_obs_max}\n", + " '''\n", + "\n", + " visits = client.query(visits_query).to_pandas().sort_values(by='obs_start_mjd')\n", + "\n", + " return visits" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "02c4c18b-fa58-4dff-8fed-a7c6030b058b", + "metadata": {}, + "outputs": [], + "source": [ + "#day_obs_min = 20250415\n", + "#day_obs_max = 20251001\n", + "day_obs_min = 20251026\n", + "#day_obs_max = 20260127\n", + "day_obs_max = 20260421\n", + "visits = fetch(day_obs_min, day_obs_max)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dce29268-b1d1-4934-889c-a10d28cc46af", + "metadata": {}, + "outputs": [], + "source": [ + "visits" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1ffb41c9-3f8d-4290-b2ca-0603b3e181f8", + "metadata": {}, + "outputs": [], + "source": [ + "visits['day_obs'].value_counts().sort_index()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "052d14c6-f2d2-49fe-ba4b-53368784feb0", + "metadata": {}, + "outputs": [], + "source": [ + "#prelsst_programs = ['BLOCK-365']\n", + "prelsst_programs = ['BLOCK-407', 'BLOCK-408', 'BLOCK-416', 'BLOCK-417', 'BLOCK-419', 'BLOCK-421']\n", + "#prelsst_programs = ['BLOCK-407']\n", + "selection_prelsst = np.in1d(visits['science_program'], prelsst_programs) \n", + "#& (visits['day_obs'] >= 20260410) * np.in1d(visits['target_name'], ['lowdust, ddf_cosmos', 'ddf_cosmos, lowdust'])\n", + "#selection_prelsst = np.in1d(visits['science_program'], prelsst_programs) & (visits['day_obs'] >= 20260410)\n", + "print(np.sum(selection_prelsst))\n", + "print(visits['day_obs'][selection_prelsst].value_counts().sort_index().to_string())\n", + "print(np.max(visits['day_obs'][selection_prelsst].value_counts()))\n", + "print(visits['day_obs'][selection_prelsst].value_counts().sort_index().to_string())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "35eaedf6-c0fc-457d-981e-6d97d4656a4c", + "metadata": {}, + "outputs": [], + "source": [ + "# 842 science program visits on 5 April\n", + "# 768 science program visits on 4 July" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b5852b2a-91d2-43c2-bf7e-49b23da17a6c", + "metadata": {}, + "outputs": [], + "source": [ + "with pd.option_context('display.max_rows', None, 'display.max_columns', None):\n", + " print(visits['day_obs'][visits['target_name'].str.contains(\"ddf_cosmos\")].value_counts().sort_index())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "aa0d9396-7c89-4422-9dc7-c7b3bef9b790", + "metadata": {}, + "outputs": [], + "source": [ + "np.unique(visits['img_type'])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2cd935f1-3f96-4cda-87b4-1edc7c9cca90", + "metadata": {}, + "outputs": [], + "source": [ + "np.sort(visits.columns)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2a241a7b-0f47-417b-9bc5-690edb119de6", + "metadata": {}, + "outputs": [], + "source": [ + "visits[selection_day_obs]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a49236c4-45c1-4e69-a7a1-dc7a9232aab1", + "metadata": {}, + "outputs": [], + "source": [ + "from astropy.time import Time" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "32192c98-f659-47dc-8722-dd683fdb45f5", + "metadata": {}, + "outputs": [], + "source": [ + "def dayObsToTime(day_obs):\n", + " year = str(day_obs)[0:4]\n", + " month = str(day_obs)[4:6]\n", + " day = str(day_obs)[6:8]\n", + " time = Time('%s-%s-%s 00:00:00'%(year, month, day), format='iso')\n", + " return time" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0d2ad00f-9999-437f-8499-b3761b991c71", + "metadata": {}, + "outputs": [], + "source": [ + "def dayObsToTwilight(day_obs, observer):\n", + " year = str(day_obs)[0:4]\n", + " month = str(day_obs)[4:6]\n", + " day = str(day_obs)[6:8]\n", + " time = Time('%s-%s-%s 23:59:00'%(year, month, day), scale='utc', format='iso')\n", + " \n", + " time_midnight = observer.midnight(time, which='nearest')\n", + " time_evening_astronomical_twilight = observer.twilight_evening_astronomical(time_midnight, which='previous')\n", + " time_morning_astronomical_twilight = observer.twilight_morning_astronomical(time_midnight, which='next')\n", + " time_evening_nautical_twilight = observer.twilight_evening_nautical(time_midnight, which='previous')\n", + " time_morning_nautical_twilight = observer.twilight_morning_nautical(time_midnight, which='next')\n", + " return time_evening_astronomical_twilight, time_morning_astronomical_twilight, time_evening_nautical_twilight, time_morning_nautical_twilight" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4daa364d-b5e1-4384-9650-3ed6ab3ef8b5", + "metadata": {}, + "outputs": [], + "source": [ + "from astroplan import Observer\n", + "observer = Observer.at_site('LSST')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fe703542-5b01-400e-9f00-311595e443fb", + "metadata": {}, + "outputs": [], + "source": [ + "import astropy.units as u" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6e29c68e-0f46-40a4-b50e-1c6c090f431c", + "metadata": {}, + "outputs": [], + "source": [ + "def convertTime(time_input, time_reference):\n", + " time_output = (time_input - (time_reference + (1. * u.day))) * u.day.to(u.hr)\n", + " return time_output.value" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dd4c5d4c-dc1a-4b34-ad9b-a750d15699c4", + "metadata": {}, + "outputs": [], + "source": [ + "\"\"\"\n", + "def gapTime(time_exp_midpt_mjd, time_evening_astronomical_twilight, time_morning_astronomical_twilight, time_day_obs):\n", + " selection_night = (time_exp_midpt_mjd > convertTime(time_evening_astronomical_twilight, time_day_obs)) & (time_exp_midpt_mjd < convertTime(time_morning_astronomical_twilight, time_day_obs))\n", + "\n", + " time_exp_midpt_mjd_augmented = np.insert(time_exp_midpt_mjd[selection_night], 0, convertTime(time_evening_astronomical_twilight, time_day_obs))\n", + " time_exp_midpt_mjd_augmented = np.append(time_exp_midpt_mjd_augmented, convertTime(time_morning_astronomical_twilight, time_day_obs))\n", + "\n", + " assert(len(time_exp_midpt_mjd[selection_night]) + 2 == len(time_exp_midpt_mjd_augmented))\n", + " assert(np.all(np.diff(time_exp_midpt_mjd_augmented) >= 0))\n", + "\n", + " selection_gap = np.diff(time_exp_midpt_mjd_augmented) * u.hr.to(u.min) > 5.\n", + "\n", + " return np.sum(np.diff(time_exp_midpt_mjd_augmented)[selection_gap])\n", + "\"\"\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4593c36b-aadc-48fc-be53-075af9f7494f", + "metadata": {}, + "outputs": [], + "source": [ + "def gapTime(time_obs_start_mjd, \n", + " time_obs_end_mjd,\n", + " time_evening_twilight, \n", + " time_morning_twilight, \n", + " time_day_obs):\n", + " \n", + " selection_night = (time_obs_end_mjd > convertTime(time_evening_twilight, time_day_obs)) & (time_obs_start_mjd < convertTime(time_morning_twilight, time_day_obs))\n", + "\n", + " gaps = time_obs_start_mjd[selection_night][1:] - time_obs_end_mjd[selection_night][0:-1]\n", + "\n", + " selection_gaps = (gaps * u.hr.to(u.min)) > 5.\n", + "\n", + " if np.any(selection_night):\n", + " gap_start_of_night = time_obs_start_mjd[selection_night][0] - convertTime(time_evening_twilight, time_day_obs)\n", + " gap_end_of_night = convertTime(time_morning_twilight, time_day_obs) - time_obs_end_mjd[selection_night][-1]\n", + " else:\n", + " # Whole night is lost\n", + " gap_start_of_night = convertTime(time_evening_twilight, time_day_obs) - convertTime(time_morning_twilight, time_day_obs)\n", + " gap_end_of_night = 0.\n", + " \n", + " return gap_start_of_night + np.sum(gaps[selection_gaps]) + gap_end_of_night" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d79dbb69-86b3-4843-9415-ae4d0c7de712", + "metadata": {}, + "outputs": [], + "source": [ + "def efficiencyFBS(time_obs_start_mjd, time_obs_end_mjd, science_program, fbs_programs):\n", + "\n", + " selection_fbs = np.isin(science_program, fbs_programs)\n", + "\n", + " gaps = time_obs_start_mjd[selection_fbs][1:] - time_obs_end_mjd[selection_fbs][0:-1]\n", + " selection_gaps = (gaps * u.hr.to(u.min)) > 5.\n", + "\n", + " if np.sum(selection_fbs) >= 2:\n", + " total_time_fbs = time_obs_end_mjd[selection_fbs][-1] - time_obs_start_mjd[selection_fbs][0]\n", + " efficiency_fbs = 1. - (np.sum(gaps[selection_gaps]) / total_time_fbs)\n", + " else:\n", + " efficiency_fbs = 0.\n", + " \n", + " return efficiency_fbs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f9275d6e-bab1-4b9a-8664-7da0b84ea9e1", + "metadata": {}, + "outputs": [], + "source": [ + "time_exp_midpt_mjd = Time(visits['exp_midpt_mjd'], format='mjd')\n", + "time_obs_start_mjd = Time(visits['obs_start_mjd'], format='mjd')\n", + "time_obs_end_mjd = Time(visits['obs_end_mjd'], format='mjd')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bbd59510-f763-4902-8da8-b74b643429a0", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(figsize=(16,8), dpi=200)\n", + "\n", + "xticks = [dayObsToTime(day_obs).mjd for day_obs in np.unique(visits['day_obs'])]\n", + "xtick_labels = np.unique(visits['day_obs'])\n", + "\n", + "day_obs_summary = {\n", + " 'day_obs': xtick_labels,\n", + " 'day_obs_mjd': xticks,\n", + " 'gap_time_nautical': [],\n", + " 'night_time_nautical': [],\n", + " 'gap_time_astronomical': [],\n", + " 'night_time_astronomical': [],\n", + " 'efficiency_fbs': [],\n", + "}\n", + "\n", + "for day_obs in np.unique(visits['day_obs']):\n", + "\n", + " time_day_obs = dayObsToTime(day_obs)\n", + " time_evening_astronomical_twilight, time_morning_astronomical_twilight, time_evening_nautical_twilight, time_morning_nautical_twilight = dayObsToTwilight(day_obs, observer)\n", + "\n", + " selection_day_obs = (visits['day_obs'] == day_obs)\n", + " #fbs_programs = ['BLOCK-365', 'BLOCK-407', 'BLOCK-408', 'BLOCK-416', 'BLOCK-417', 'BLOCK-419', 'BLOCK-421'] # Commissioning\n", + " fbs_programs = ['BLOCK-407', 'BLOCK-408', 'BLOCK-416', 'BLOCK-417', 'BLOCK-419', 'BLOCK-421'] # Early Ops\n", + " aos_programs = ['BLOCK-T539']\n", + "\n", + " time_exp_midpt_mjd = convertTime(\n", + " Time(visits['exp_midpt_mjd'][selection_day_obs], format='mjd'), \n", + " time_day_obs\n", + " )\n", + " time_obs_start_mjd = convertTime(\n", + " Time(visits['exp_midpt_mjd'][selection_day_obs], format='mjd'), \n", + " time_day_obs\n", + " )\n", + " time_obs_end_mjd = convertTime(\n", + " Time(visits['exp_midpt_mjd'][selection_day_obs], format='mjd'), \n", + " time_day_obs\n", + " )\n", + " x = np.tile(time_day_obs.mjd, len(visits[selection_day_obs]))\n", + "\n", + " selection_fbs = np.isin(visits[selection_day_obs]['science_program'], fbs_programs)\n", + " selection_aos = np.isin(visits[selection_day_obs]['science_program'], aos_programs)\n", + "\n", + " plt.scatter(x[~selection_fbs], time_exp_midpt_mjd[~selection_fbs], marker='_', c='tab:orange')\n", + " plt.scatter(x[selection_fbs], time_exp_midpt_mjd[selection_fbs], marker='_', c='tab:blue')\n", + " plt.scatter(x[selection_aos], time_exp_midpt_mjd[selection_aos], marker='_', c='tab:green')\n", + "\n", + " plt.scatter(time_day_obs.mjd, convertTime(time_evening_astronomical_twilight, time_day_obs), marker='o', c='black')\n", + " plt.scatter(time_day_obs.mjd, convertTime(time_morning_astronomical_twilight, time_day_obs), marker='o', c='black')\n", + " plt.scatter(time_day_obs.mjd, convertTime(time_evening_nautical_twilight, time_day_obs), marker='o', c='0.5')\n", + " plt.scatter(time_day_obs.mjd, convertTime(time_morning_nautical_twilight, time_day_obs), marker='o', c='0.5')\n", + " #plt.plot(time_day_obs.mjd, convertTime(time_evening_astronomical_twilight, time_day_obs))\n", + "\n", + " # Astronomical twilight\n", + " #day_obs_summary['gap_time'].append(gapTime(time_exp_midpt_mjd, time_evening_astronomical_twilight, time_morning_astronomical_twilight, time_day_obs))\n", + " #day_obs_summary['night_time'].append((time_morning_astronomical_twilight - time_evening_astronomical_twilight).value * u.day.to(u.hr))\n", + "\n", + " # Nautical twilight\n", + " #day_obs_summary['gap_time'].append(gapTime(time_exp_midpt_mjd, time_evening_nautical_twilight, time_morning_nautical_twilight, time_day_obs))\n", + " #day_obs_summary['night_time'].append((time_morning_nautical_twilight - time_evening_nautical_twilight).value * u.day.to(u.hr))\n", + "\n", + " day_obs_summary['gap_time_nautical'].append(gapTime(time_obs_start_mjd, time_obs_end_mjd, \n", + " time_evening_nautical_twilight, time_morning_nautical_twilight, \n", + " time_day_obs))\n", + " day_obs_summary['night_time_nautical'].append((time_morning_nautical_twilight - time_evening_nautical_twilight).value * u.day.to(u.hr))\n", + " day_obs_summary['gap_time_astronomical'].append(gapTime(time_obs_start_mjd, time_obs_end_mjd, \n", + " time_evening_astronomical_twilight, time_morning_astronomical_twilight, \n", + " time_day_obs))\n", + " day_obs_summary['night_time_astronomical'].append((time_morning_astronomical_twilight - time_evening_astronomical_twilight).value * u.day.to(u.hr))\n", + "\n", + " day_obs_summary['efficiency_fbs'].append(efficiencyFBS(time_obs_start_mjd, time_obs_end_mjd, visits[selection_day_obs]['science_program'], fbs_programs))\n", + "\n", + "plt.ylim(-1., 10.5) # Early Ops\n", + "#plt.ylim(-2., 11.) # Commissioning\n", + "plt.xticks(xticks, xtick_labels, rotation=90., fontsize=8)\n", + "plt.ylabel('Time (hrs)')\n", + "plt.tight_layout()\n", + "plt.show()\n", + "\n", + "for key in day_obs_summary.keys():\n", + " day_obs_summary[key] = np.array([_ for _ in day_obs_summary[key]])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3e096a9a-42d1-4d31-a689-fce6fea1f2dd", + "metadata": {}, + "outputs": [], + "source": [ + "#gapTime(time_exp_midpt_mjd, time_evening_astronomical_twilight, time_morning_astronomical_twilight, time_day_obs)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2f3e8513-21a9-42e5-89e6-06feb1cf0392", + "metadata": {}, + "outputs": [], + "source": [ + "time_evening_astronomical_twilight.mjd" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b752c663-5230-460a-9368-414f16478e6e", + "metadata": {}, + "outputs": [], + "source": [ + "(time_morning_astronomical_twilight - time_evening_astronomical_twilight).value * u.day.to(u.hr)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2e9e30fa-850f-4845-82f9-5deac21d0449", + "metadata": {}, + "outputs": [], + "source": [ + "day_obs_summary" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bc61be40-88af-4d40-bc9a-9d7924ac88c8", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(figsize=(12,8))\n", + "plt.scatter(day_obs_summary['day_obs_mjd'], 1. - (day_obs_summary['gap_time_astronomical'] / day_obs_summary['night_time_astronomical']), label='Between Astronomical Twilights')\n", + "plt.scatter(day_obs_summary['day_obs_mjd'], 1. - (day_obs_summary['gap_time_nautical'] / day_obs_summary['night_time_nautical']), label='Between Nautical Twilights')\n", + "plt.scatter(day_obs_summary['day_obs_mjd'], day_obs_summary['efficiency_fbs'], label='Within FBS pre-LSST')\n", + "plt.grid()\n", + "plt.ylim(0., 1.)\n", + "plt.xticks(xticks, xtick_labels, rotation=90.)\n", + "#plt.ylabel('Fraction Inactive\\n(sum of visit gaps > 5 min)')\n", + "plt.ylabel('Active Observing Fraction\\n(visit gaps < 5 min)')\n", + "plt.legend()\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fed95603-6da2-4921-823f-e593127f55a9", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "plt.hist(day_obs_summary['gap_time'] / day_obs_summary['night_time'], bins=np.linspace(0., 1., 21))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a4426ccd-f34e-4f51-97d2-a28b18425e9d", + "metadata": {}, + "outputs": [], + "source": [ + "np.mean(day_obs_summary['gap_time'] / day_obs_summary['night_time'])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8f162ec1-fad8-4f38-94a7-9b4362e4e3c1", + "metadata": {}, + "outputs": [], + "source": [ + "for _ in " + ] + }, + { + "cell_type": "markdown", + "id": "e7583b94-77a7-4a37-a6b6-df060d386fe6", + "metadata": {}, + "source": [ + "# Visit Gap Time Distribution" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4f16e98e-f7df-4455-a806-6a55fc2d0989", + "metadata": {}, + "outputs": [], + "source": [ + "fbs_programs = ['BLOCK-407', 'BLOCK-408', 'BLOCK-416', 'BLOCK-417', 'BLOCK-419']\n", + "selection_fbs = np.isin(visits['science_program'], fbs_programs)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c772e681-f2eb-428d-92a1-bccb6e72a323", + "metadata": {}, + "outputs": [], + "source": [ + "visit_gaps = pd.DataFrame({\n", + " 'day_obs': visits[selection_fbs]['day_obs'][1:],\n", + " 'visit_gap': (visits[selection_fbs]['obs_start_mjd'][1:].values - visits[selection_fbs]['obs_end_mjd'][0:-1].values) * u.day.to(u.second)\n", + "})\n", + "\n", + "\n", + "groups = visit_gaps.groupby('day_obs')\n", + "visit_gaps_summary = pd.DataFrame({\n", + " 'day_obs': groups['day_obs'].first(),\n", + " 'n_visits': groups['day_obs'].count(),\n", + " 'visit_gap_10': groups['visit_gap'].quantile(0.1),\n", + " 'visit_gap_50': groups['visit_gap'].quantile(0.50),\n", + " 'visit_gap_90': groups['visit_gap'].quantile(0.9),\n", + "})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "56c4a36c-73d0-477d-a90e-a36d04d5bd85", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "plt.scatter(visit_gaps_summary['day_obs'], visit_gaps_summary['visit_gap_50'])\n", + "plt.ylim(0., 30.)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9cf5fd6f-8811-46be-813f-5f873977fe10", + "metadata": {}, + "outputs": [], + "source": [ + "bins = np.linspace(0., 30., 31)\n", + "\n", + "#day_obs_threshold = 20260109\n", + "day_obs_threshold = 20260401\n", + "\n", + "#selection_day_obs = (visit_gaps['day_obs'] >= 20251203)\n", + "selection_day_obs = (visit_gaps['day_obs'] >= day_obs_threshold)\n", + "\n", + "\n", + "\n", + "plt.figure()\n", + "#plt.hist(visit_gaps['visit_gap'][~selection_day_obs], bins=bins, histtype='step', density=True, label=f\"day_obs < {day_obs_threshold}\")\n", + "plt.hist(visit_gaps['visit_gap'][selection_day_obs], bins=bins, histtype='step', density=True, label=f\"day_obs >= {day_obs_threshold}\")\n", + "\n", + "#plt.axvline(np.median(visit_gaps['visit_gap'][~selection_day_obs]))\n", + "#plt.axvline(np.median(visit_gaps['visit_gap'][selection_day_obs]))\n", + "\n", + "#plt.xlim(0., 60.)\n", + "plt.xlim(0., 30.)\n", + "plt.legend()\n", + "plt.xlabel('Visit Gap (second)')\n", + "plt.ylabel('PDF')\n", + "plt.title('FBS pre-LSST Visits')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7291cf2c-1f3e-4b8f-bf23-9bf05098e3f1", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "LSST", + "language": "python", + "name": "lsst" + }, + "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.13.9" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/notebooks/image_quality_minimal.ipynb b/notebooks/image_quality_minimal.ipynb new file mode 100644 index 0000000..040acb6 --- /dev/null +++ b/notebooks/image_quality_minimal.ipynb @@ -0,0 +1,299 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "59992f52-3886-44d7-b947-0ea83f74b26b", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import pandas as pd\n", + "\n", + "import matplotlib.pyplot as plt\n", + "\n", + "%matplotlib widget" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3a3fd2c1-4141-4de3-a1d8-dd2563bd7d91", + "metadata": {}, + "outputs": [], + "source": [ + "from lsst.utils.plotting import publication_plots\n", + "from lsst.utils.plotting import get_multiband_plot_colors\n", + "\n", + "colors = get_multiband_plot_colors()\n", + "bands = colors.keys()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1f9ed8ef-6652-45f8-a8dd-9e017f39d4e3", + "metadata": {}, + "outputs": [], + "source": [ + "from lsst.summit.utils import (\n", + " getAirmassSeeingCorrection,\n", + " getBandpassSeeingCorrection,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fc581d34-90bb-462e-ac2d-bbfbdff2cf87", + "metadata": {}, + "outputs": [], + "source": [ + "CAM_FWHM = 0.207 # arcsec" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "177bee50-c745-4d16-979f-2b2fa38a3971", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "os.environ[\"no_proxy\"] += \",.consdb\"\n", + "\n", + "from lsst.summit.utils import ConsDbClient\n", + "\n", + "client = ConsDbClient(\"http://consdb-pq.consdb:8080/consdb\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bba9fded-2295-4534-b8c7-e73bbc422321", + "metadata": {}, + "outputs": [], + "source": [ + "instrument = 'lsstcam'" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "47f5cc04-52f3-4932-b8b4-984683a0dda0", + "metadata": {}, + "outputs": [], + "source": [ + "def fetch(day_obs_min, day_obs_max):\n", + " \n", + " visits_query = f'''\n", + " SELECT \n", + " ccdvisit1_quicklook.psf_area,\n", + " ccdvisit1_quicklook.psf_sigma,\n", + " ccdvisit1_quicklook.psf_ixx,\n", + " ccdvisit1_quicklook.psf_ixy,\n", + " ccdvisit1_quicklook.psf_iyy,\n", + " ccdvisit1.detector,\n", + " visit1.visit_id,\n", + " visit1.seq_num,\n", + " visit1.band,\n", + " visit1.physical_filter,\n", + " visit1.day_obs,\n", + " visit1.target_name,\n", + " visit1.science_program,\n", + " visit1.observation_reason,\n", + " visit1.airmass,\n", + " visit1.altitude,\n", + " visit1.azimuth,\n", + " visit1.sky_rotation,\n", + " visit1.exp_midpt_mjd,\n", + " visit1.dimm_seeing,\n", + " visit1.s_ra,\n", + " visit1.s_dec,\n", + " visit1_quicklook.psf_sigma_min,\n", + " visit1_quicklook.psf_sigma_median,\n", + " visit1_quicklook.aos_fwhm,\n", + " visit1_quicklook.donut_blur_fwhm,\n", + " visit1_quicklook.ringss_seeing,\n", + " visit1_quicklook.seeing_zenith_500nm_min,\n", + " visit1_quicklook.seeing_zenith_500nm_median,\n", + " visit1_quicklook.physical_rotator_angle\n", + " FROM\n", + " cdb_{instrument}.ccdvisit1_quicklook AS ccdvisit1_quicklook,\n", + " cdb_{instrument}.ccdvisit1 AS ccdvisit1,\n", + " cdb_{instrument}.visit1_quicklook AS visit1_quicklook,\n", + " cdb_{instrument}.visit1 AS visit1 \n", + " WHERE \n", + " ccdvisit1.ccdvisit_id = ccdvisit1_quicklook.ccdvisit_id\n", + " AND ccdvisit1.visit_id = visit1.visit_id \n", + " AND visit1.visit_id = visit1_quicklook.visit_id\n", + " AND ccdvisit1.detector NOT IN (168, 188, 123, 27, 0, 20, 65, 161) -- vignetted\n", + " AND ccdvisit1.detector NOT IN (191, 192, 195, 196, 199, 200, 203, 204) -- corner wavefront sensors\n", + " AND visit1.airmass > 0\n", + " AND visit1.day_obs >= {day_obs_min} AND visit1.day_obs <= {day_obs_max}\n", + " AND visit1.science_program in ('BLOCK-365', 'BLOCK-407', 'BLOCK-408', 'BLOCK-416', 'BLOCK-417', 'BLOCK-419');\n", + " '''\n", + " \n", + " ccdvisits = client.query(visits_query).to_pandas()\n", + "\n", + " pixel_scale = 0.2\n", + " sig2fwhm = 2 * np.sqrt(2 * np.log(2))\n", + " ccdvisits[\"psf_fwhm\"] = ccdvisits[\"psf_sigma\"] * sig2fwhm * pixel_scale\n", + " ccdvisits[\"psf_fwhm\"] = pd.to_numeric(ccdvisits[\"psf_fwhm\"], errors=\"coerce\")\n", + " ccdvisits[\"donut_blur_fwhm\"] = pd.to_numeric(ccdvisits[\"donut_blur_fwhm\"], errors=\"coerce\")\n", + " ccdvisits[\"aos_fwhm\"] = pd.to_numeric(ccdvisits[\"aos_fwhm\"], errors=\"coerce\")\n", + "\n", + " ccdvisits[\"psf_e1\"] = (ccdvisits[\"psf_ixx\"] - ccdvisits[\"psf_iyy\"]) / (ccdvisits[\"psf_ixx\"] + ccdvisits[\"psf_iyy\"])\n", + " ccdvisits[\"psf_e2\"] = (2. * ccdvisits[\"psf_ixy\"]) / (ccdvisits[\"psf_ixx\"] + ccdvisits[\"psf_iyy\"])\n", + " ccdvisits[\"psf_e\"] = np.sqrt(np.array(ccdvisits[\"psf_e1\"]**2 + ccdvisits[\"psf_e2\"]**2, dtype=float))\n", + "\n", + " ccdvisits[\"psf_e\"] = pd.to_numeric(ccdvisits[\"psf_e\"], errors=\"coerce\")\n", + " ccdvisits[\"psf_e1\"] = pd.to_numeric(ccdvisits[\"psf_e1\"], errors=\"coerce\")\n", + " ccdvisits[\"psf_e2\"] = pd.to_numeric(ccdvisits[\"psf_e2\"], errors=\"coerce\")\n", + "\n", + " return ccdvisits" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "041bcfb4-c28d-41ae-9dff-20bea3a8b483", + "metadata": {}, + "outputs": [], + "source": [ + "day_obs_min = 20250119\n", + "day_obs_max = 20260119\n", + "ccdvisits = fetch(day_obs_min, day_obs_max)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b5714f29-5ddf-4f2f-8d9d-21ae4dfc7b00", + "metadata": {}, + "outputs": [], + "source": [ + "airmass_seeing_correction = np.array([getAirmassSeeingCorrection(airmass) for airmass in ccdvisits[\"airmass\"]])\n", + "bandpass_seeing_correction = np.array([getBandpassSeeingCorrection(band) for band in ccdvisits[\"physical_filter\"]])\n", + "ccdvisits[\"psf_fwhm_zenith_500nm\"] = ccdvisits[\"psf_fwhm\"] * airmass_seeing_correction * bandpass_seeing_correction" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8be116b3-cdb7-4137-9eec-1c0411cb9582", + "metadata": {}, + "outputs": [], + "source": [ + "ccdvisits['donut_blur_atm_fwhm'] = np.sqrt(ccdvisits['donut_blur_fwhm']**2 - CAM_FWHM**2)\n", + "ccdvisits['aos_cam_fwhm'] = np.sqrt(ccdvisits['aos_fwhm']**2 + CAM_FWHM**2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9ccffd97-0014-4f5c-a08f-3c474dd26934", + "metadata": {}, + "outputs": [], + "source": [ + "groups = ccdvisits.groupby('visit_id')\n", + "visits_summary = pd.DataFrame({\n", + " 'day_obs': groups['day_obs'].first(),\n", + " 'target_name': groups['target_name'].first(),\n", + " 'science_program': groups['science_program'].first(),\n", + " 'observation_reason': groups['observation_reason'].first(),\n", + " 'seq_num': groups['seq_num'].median(),\n", + " 'exp_midpt_mjd': groups['exp_midpt_mjd'].median(),\n", + " 'donut_blur_fwhm': groups['donut_blur_fwhm'].median(),\n", + " 'aos_fwhm': groups['aos_fwhm'].median(),\n", + " 'donut_blur_atm_fwhm': groups['donut_blur_atm_fwhm'].median(),\n", + " 'aos_cam_fwhm': groups['aos_cam_fwhm'].median(),\n", + " 'physical_rotator_angle': groups['physical_rotator_angle'].median(),\n", + " 'altitude': groups['altitude'].median(),\n", + " 'psf_fwhm_05': groups['psf_fwhm'].quantile(0.05),\n", + " 'psf_fwhm_50': groups['psf_fwhm'].quantile(0.50),\n", + " 'psf_fwhm_95': groups['psf_fwhm'].quantile(0.95),\n", + " 'psf_fwhm_zenith_500nm_50': groups['psf_fwhm_zenith_500nm'].quantile(0.50),\n", + " 'psf_e_05': groups['psf_e'].quantile(0.05),\n", + " 'psf_e_50': groups['psf_e'].quantile(0.50),\n", + " 'psf_e_95': groups['psf_e'].quantile(0.95),\n", + " 'psf_e1_05': groups['psf_e1'].quantile(0.05),\n", + " 'psf_e1_50': groups['psf_e1'].quantile(0.50),\n", + " 'psf_e1_95': groups['psf_e1'].quantile(0.95),\n", + " 'psf_e2_05': groups['psf_e2'].quantile(0.05),\n", + " 'psf_e2_50': groups['psf_e2'].quantile(0.50),\n", + " 'psf_e2_95': groups['psf_e2'].quantile(0.95),\n", + " 'band': groups['band'].first(),\n", + "})\n", + "visits_summary['psf_fwhm_95_05'] = np.sqrt(visits_summary['psf_fwhm_95']**2 - visits_summary['psf_fwhm_05']**2)\n", + "visits_summary['sys_50'] = np.sqrt(visits_summary['psf_fwhm_50']**2 + CAM_FWHM**2 - visits_summary['donut_blur_fwhm']**2)\n", + "visits_summary['psf_fwhm_model'] = np.sqrt(visits_summary['aos_fwhm']**2 + visits_summary['donut_blur_fwhm']**2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21429500-7988-4e9a-b8ec-61b834d519db", + "metadata": {}, + "outputs": [], + "source": [ + "# Exclude y due to issues with estimating donut blur\n", + "selection = (visits_summary['band'] != \"y\") & ~visits_summary['observation_reason'].str.contains('filter_change_close_loop')\n", + "\n", + "groups = visits_summary[selection].groupby('day_obs')\n", + "day_obs_summary = pd.DataFrame({\n", + " 'day_obs': groups['day_obs'].first(),\n", + " 'n_visits': groups['day_obs'].count(),\n", + " 'psf_fwhm_95_05_low': groups['psf_fwhm_95_05'].quantile(0.10),\n", + " 'psf_fwhm_95_05_50': groups['psf_fwhm_95_05'].quantile(0.50),\n", + " 'psf_fwhm_95_05_high': groups['psf_fwhm_95_05'].quantile(0.90),\n", + " 'aos_fwhm_low': groups['aos_fwhm'].quantile(0.10),\n", + " 'aos_fwhm_50': groups['aos_fwhm'].quantile(0.50),\n", + " 'aos_fwhm_high': groups['aos_fwhm'].quantile(0.90),\n", + " 'aos_cam_fwhm_low': groups['aos_cam_fwhm'].quantile(0.10),\n", + " 'aos_cam_fwhm_50': groups['aos_cam_fwhm'].quantile(0.50),\n", + " 'aos_cam_fwhm_high': groups['aos_cam_fwhm'].quantile(0.90),\n", + " 'sys_50_low': groups['sys_50'].quantile(0.10),\n", + " 'sys_50_50': groups['sys_50'].quantile(0.50),\n", + " 'sys_50_high': groups['sys_50'].quantile(0.90),\n", + " 'psf_e_50_low': groups['psf_e_50'].quantile(0.10),\n", + " 'psf_e_50_50': groups['psf_e_50'].quantile(0.50),\n", + " 'psf_e_50_high': groups['psf_e_50'].quantile(0.90),\n", + " 'psf_e1_50_low': groups['psf_e1_50'].quantile(0.10),\n", + " 'psf_e1_50_50': groups['psf_e1_50'].quantile(0.50),\n", + " 'psf_e1_50_high': groups['psf_e1_50'].quantile(0.90),\n", + " 'psf_e2_50_low': groups['psf_e2_50'].quantile(0.10),\n", + " 'psf_e2_50_50': groups['psf_e2_50'].quantile(0.50),\n", + " 'psf_e2_50_high': groups['psf_e2_50'].quantile(0.90),\n", + " 'psf_fwhm_50_low': groups['psf_fwhm_50'].quantile(0.10),\n", + " 'psf_fwhm_50_50': groups['psf_fwhm_50'].quantile(0.50),\n", + " 'psf_fwhm_50_high': groups['psf_fwhm_50'].quantile(0.90),\n", + " 'donut_blur_atm_fwhm_low': groups['donut_blur_atm_fwhm'].quantile(0.10),\n", + " 'donut_blur_atm_fwhm_50': groups['donut_blur_atm_fwhm'].quantile(0.50),\n", + " 'donut_blur_atm_fwhm_high': groups['donut_blur_atm_fwhm'].quantile(0.90),\n", + "})" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "LSST", + "language": "python", + "name": "lsst" + }, + "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.13.9" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/notebooks/image_quality_trending.ipynb b/notebooks/image_quality_trending.ipynb new file mode 100644 index 0000000..a8c4b84 --- /dev/null +++ b/notebooks/image_quality_trending.ipynb @@ -0,0 +1,4997 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "11e90e45-38f8-44f5-93df-e0a0b81b73af", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import pandas as pd\n", + "\n", + "import matplotlib.pyplot as plt\n", + "\n", + "#%matplotlib widget" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8a6fd5a8-9017-41ea-9f45-17a399657e77", + "metadata": {}, + "outputs": [], + "source": [ + "from lsst.utils.plotting import publication_plots\n", + "from lsst.utils.plotting import get_multiband_plot_colors\n", + "\n", + "colors = get_multiband_plot_colors()\n", + "bands = colors.keys()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7150d061-b664-4007-a10a-c88146dbd21f", + "metadata": {}, + "outputs": [], + "source": [ + "from lsst.summit.utils import (\n", + " getAirmassSeeingCorrection,\n", + " getBandpassSeeingCorrection,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b47e74a1-645f-41eb-896f-399f1c762d19", + "metadata": {}, + "outputs": [], + "source": [ + "CAM_FWHM = 0.207 # arcsec" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2add4d51-afec-4114-b053-76ac9024aaa4", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "os.environ[\"no_proxy\"] += \",.consdb\"\n", + "\n", + "from lsst.summit.utils import ConsDbClient\n", + "\n", + "client = ConsDbClient(\"http://consdb-pq.consdb:8080/consdb\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a53dc0f4-4094-4322-992c-0638574d2048", + "metadata": {}, + "outputs": [], + "source": [ + "instrument = 'lsstcam'" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5dd1e4fb-25d8-4cac-aa17-b77156e652ec", + "metadata": {}, + "outputs": [], + "source": [ + "def fetch(day_obs_min, day_obs_max):\n", + " \n", + " visits_query = f'''\n", + " SELECT \n", + " ccdvisit1_quicklook.psf_area,\n", + " ccdvisit1_quicklook.psf_sigma,\n", + " ccdvisit1_quicklook.psf_ixx,\n", + " ccdvisit1_quicklook.psf_ixy,\n", + " ccdvisit1_quicklook.psf_iyy,\n", + " --ccdvisit1_quicklook.z4,\n", + " --ccdvisit1_quicklook.z5,\n", + " --ccdvisit1_quicklook.z6,\n", + " --ccdvisit1_quicklook.z7,\n", + " --ccdvisit1_quicklook.z8,\n", + " --ccdvisit1_quicklook.z9,\n", + " --ccdvisit1_quicklook.z10,\n", + " --ccdvisit1_quicklook.z11,\n", + " --ccdvisit1_quicklook.z12,\n", + " --ccdvisit1_quicklook.z13,\n", + " --ccdvisit1_quicklook.z14,\n", + " --ccdvisit1_quicklook.z15,\n", + " --ccdvisit1_quicklook.z16,\n", + " --ccdvisit1_quicklook.z17,\n", + " --ccdvisit1_quicklook.z18,\n", + " --ccdvisit1_quicklook.z19,\n", + " --ccdvisit1_quicklook.z20,\n", + " --ccdvisit1_quicklook.z21,\n", + " --ccdvisit1_quicklook.z22,\n", + " --ccdvisit1_quicklook.z23,\n", + " --ccdvisit1_quicklook.z24,\n", + " --ccdvisit1_quicklook.z25,\n", + " --ccdvisit1_quicklook.z26,\n", + " --ccdvisit1_quicklook.z27,\n", + " --ccdvisit1_quicklook.z28,\n", + " ccdvisit1.detector,\n", + " visit1.visit_id,\n", + " visit1.seq_num,\n", + " visit1.band,\n", + " visit1.physical_filter,\n", + " visit1.day_obs,\n", + " visit1.target_name,\n", + " visit1.science_program,\n", + " visit1.observation_reason,\n", + " visit1.airmass,\n", + " visit1.altitude,\n", + " visit1.azimuth,\n", + " visit1.sky_rotation,\n", + " visit1.exp_midpt_mjd,\n", + " visit1.dimm_seeing,\n", + " visit1.s_ra,\n", + " visit1.s_dec,\n", + " visit1_quicklook.psf_sigma_min,\n", + " visit1_quicklook.psf_sigma_median,\n", + " visit1_quicklook.aos_fwhm,\n", + " visit1_quicklook.donut_blur_fwhm,\n", + " visit1_quicklook.ringss_seeing,\n", + " visit1_quicklook.seeing_zenith_500nm_min,\n", + " visit1_quicklook.seeing_zenith_500nm_median,\n", + " visit1_quicklook.physical_rotator_angle\n", + " FROM\n", + " cdb_{instrument}.ccdvisit1_quicklook AS ccdvisit1_quicklook,\n", + " cdb_{instrument}.ccdvisit1 AS ccdvisit1,\n", + " cdb_{instrument}.visit1_quicklook AS visit1_quicklook,\n", + " cdb_{instrument}.visit1 AS visit1 \n", + " WHERE \n", + " ccdvisit1.ccdvisit_id = ccdvisit1_quicklook.ccdvisit_id\n", + " AND ccdvisit1.visit_id = visit1.visit_id \n", + " AND visit1.visit_id = visit1_quicklook.visit_id\n", + " AND ccdvisit1.detector NOT IN (168, 188, 123, 27, 0, 20, 65, 161) -- vignetted\n", + " AND ccdvisit1.detector NOT IN (191, 192, 195, 196, 199, 200, 203, 204) -- corner wavefront sensors\n", + " AND visit1.airmass > 0\n", + " AND visit1.day_obs >= {day_obs_min} AND visit1.day_obs <= {day_obs_max}\n", + " -- AND visit1.science_program in ('BLOCK-365', 'BLOCK-407');\n", + " AND visit1.science_program in ('BLOCK-365', 'BLOCK-407', 'BLOCK-408', 'BLOCK-416', 'BLOCK-417', 'BLOCK-419', 'BLOCK-421') \n", + " -- , 'BLOCK-T698')\n", + " -- AND visit1.science_program in ('BLOCK-T539', 'BLOCK-T698', 'BLOCK-T641') -- initial alignment, closed loop stability\n", + " AND visit1.physical_filter in ('u_24', 'g_6', 'r_57', 'i_39', 'z_20', 'y_10');\n", + " '''\n", + " \n", + " ccdvisits = client.query(visits_query).to_pandas()\n", + "\n", + " pixel_scale = 0.2\n", + " sig2fwhm = 2 * np.sqrt(2 * np.log(2))\n", + " ccdvisits[\"psf_fwhm\"] = ccdvisits[\"psf_sigma\"] * sig2fwhm * pixel_scale\n", + " ccdvisits[\"psf_fwhm\"] = pd.to_numeric(ccdvisits[\"psf_fwhm\"], errors=\"coerce\")\n", + " ccdvisits[\"donut_blur_fwhm\"] = pd.to_numeric(ccdvisits[\"donut_blur_fwhm\"], errors=\"coerce\")\n", + " ccdvisits[\"aos_fwhm\"] = pd.to_numeric(ccdvisits[\"aos_fwhm\"], errors=\"coerce\")\n", + "\n", + " ccdvisits[\"psf_area\"] = pd.to_numeric(ccdvisits[\"psf_area\"], errors=\"coerce\")\n", + " ccdvisits[\"psf_fwhm_area\"] = 0.663 * pixel_scale * np.sqrt(ccdvisits[\"psf_area\"])\n", + "\n", + " #airmass_seeing_correction = np.array([getAirmassSeeingCorrection(airmass) for airmass in ccdvisits[\"airmass\"]])\n", + " #bandpass_seeing_correction = np.array([getBandpassSeeingCorrection(band) for band in ccdvisits[\"physical_filter\"]])\n", + " #ccdvisits[\"psf_fwhm_zenith_500nm\"] = ccdvisits[\"psf_fwhm\"] * airmass_seeing_correction * bandpass_seeing_correction\n", + "\n", + " #ccdvisits[\"psf_e1\"] = (ccdvisits[\"psf_ixx\"]**2 - ccdvisits[\"psf_iyy\"]**2) / (ccdvisits[\"psf_ixx\"]**2 + ccdvisits[\"psf_iyy\"]**2)\n", + " #ccdvisits[\"psf_e2\"] = (2. * ccdvisits[\"psf_ixy\"]**2) / (ccdvisits[\"psf_ixx\"]**2 + ccdvisits[\"psf_iyy\"]**2)\n", + " ccdvisits[\"psf_e1\"] = (ccdvisits[\"psf_ixx\"] - ccdvisits[\"psf_iyy\"]) / (ccdvisits[\"psf_ixx\"] + ccdvisits[\"psf_iyy\"])\n", + " ccdvisits[\"psf_e2\"] = (2. * ccdvisits[\"psf_ixy\"]) / (ccdvisits[\"psf_ixx\"] + ccdvisits[\"psf_iyy\"])\n", + " ccdvisits[\"psf_e\"] = np.sqrt(np.array(ccdvisits[\"psf_e1\"]**2 + ccdvisits[\"psf_e2\"]**2, dtype=float))\n", + " ccdvisits[\"psf_e1_unnormalized\"] = (ccdvisits[\"psf_ixx\"] - ccdvisits[\"psf_iyy\"])\n", + " ccdvisits[\"psf_e2_unnormalized\"] = (2. * ccdvisits[\"psf_ixy\"])\n", + " ccdvisits[\"psf_e_unnormalized\"] = np.sqrt(np.array(ccdvisits[\"psf_e1_unnormalized\"]**2 + ccdvisits[\"psf_e2_unnormalized\"]**2, dtype=float))\n", + "\n", + " ccdvisits[\"psf_e\"] = pd.to_numeric(ccdvisits[\"psf_e\"], errors=\"coerce\")\n", + " ccdvisits[\"psf_e1\"] = pd.to_numeric(ccdvisits[\"psf_e1\"], errors=\"coerce\")\n", + " ccdvisits[\"psf_e2\"] = pd.to_numeric(ccdvisits[\"psf_e2\"], errors=\"coerce\")\n", + "\n", + " return ccdvisits" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ef9d10ac-dce3-4686-b932-19d09b832b29", + "metadata": {}, + "outputs": [], + "source": [ + "#day_obs_min = 20250620 # Start of SV surveys\n", + "#day_obs_max = 20250709\n", + "#day_obs_min = 20250906\n", + "#day_obs_max = 20250906\n", + "#day_obs_min = 20250708\n", + "#day_obs_max = 20250708\n", + "#day_obs_min = 20250709\n", + "#day_obs_max = 20250709\n", + "#day_obs_min = 20250913\n", + "#day_obs_max = 20250913\n", + "#day_obs_min = 20250915\n", + "#day_obs_max = 20250915\n", + "#day_obs_min = 20250921\n", + "#day_obs_max = 20250921 # Last night of SV surveys\n", + "#day_obs_min = 20251026\n", + "#day_obs_max = 20251026\n", + "#day_obs_min = 20251027\n", + "#day_obs_max = 20251027\n", + "#day_obs_min = 20251029\n", + "#day_obs_max = 20251029\n", + "#day_obs_max = 20251102\n", + "#day_obs_min = 20251104\n", + "#day_obs_max = 20251125\n", + "#day_obs_max = 20251201\n", + "#day_obs_min = 20251113\n", + "#day_obs_min = 20251126\n", + "#day_obs_min = 20251128\n", + "#day_obs_min = 20251230 # Start of sustained pre-LSST observations\n", + "#day_obs_min = 20260106\n", + "#day_obs_min = 20251230\n", + "#day_obs_max = 20260107\n", + "#day_obs_min = 20260113\n", + "#day_obs_max = 20260127\n", + "#day_obs_min = 20260219\n", + "#day_obs_max = 20260223\n", + "#day_obs_min = 20260224\n", + "#day_obs_min = 20260308\n", + "#day_obs_max = 20260309 # End of sustained pre-LSST observations\n", + "#day_obs_min = 20260328\n", + "#day_obs_max = 20260330\n", + "day_obs_min = 20260401\n", + "#day_obs_min = 20260402\n", + "#day_obs_max = 20260407\n", + "#day_obs_min = 20260414\n", + "day_obs_max = 20260423\n", + "ccdvisits = fetch(day_obs_min, day_obs_max)\n", + "print(len(ccdvisits))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "227549af-1ca9-4af8-b46a-9b85cf379ba1", + "metadata": {}, + "outputs": [], + "source": [ + "np.unique(ccdvisits[\"physical_filter\"])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8f15686a-7301-4c2e-a225-65b54d1a2cd8", + "metadata": {}, + "outputs": [], + "source": [ + "np.unique(ccdvisits[\"science_program\"])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "caf757f9-ee21-46a6-b0fc-ff19555611f2", + "metadata": {}, + "outputs": [], + "source": [ + "airmass_seeing_correction = np.array([getAirmassSeeingCorrection(airmass) for airmass in ccdvisits[\"airmass\"]])\n", + "bandpass_seeing_correction = np.array([getBandpassSeeingCorrection(band) for band in ccdvisits[\"physical_filter\"]])\n", + "ccdvisits[\"psf_fwhm_zenith_500nm\"] = ccdvisits[\"psf_fwhm\"] * airmass_seeing_correction * bandpass_seeing_correction" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "092e3e24-5800-4a7b-b87b-18594bf8f4fd", + "metadata": {}, + "outputs": [], + "source": [ + "### IMPROVED METHOD ###\n", + "\n", + "# fwhm_geom = visits[psf_col] * SIGMA_TO_FWHM * pixel_scale\n", + "# idiq_aos_cam = np.sqrt(visits.aos_fwhm**2 + CAM_FWHM**2)\n", + "# atm_fwhm_aos_cam = np.sqrt(fwhm_geom**2 - idiq_aos_cam**2)\n", + "# atm_500_zenith = atm_fwhm_aos_cam / wavelen_corrections / airmass_corrections\n", + "# fwhm_500_zenith = np.sqrt(atm_500_zenith**2 + idiq_aos_cam**2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1b56912c-ed27-4b57-8a5f-d10a815c609b", + "metadata": {}, + "outputs": [], + "source": [ + "ccdvisits['donut_blur_atm_fwhm'] = np.sqrt(ccdvisits['donut_blur_fwhm']**2 - CAM_FWHM**2)\n", + "ccdvisits['aos_cam_fwhm'] = np.sqrt(ccdvisits['aos_fwhm']**2 + CAM_FWHM**2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ffa4926b-afa2-4f3a-b922-11469c7a239a", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(figsize=(6, 6))\n", + "#plt.scatter(ccdvisits['psf_fwhm'], ccdvisits['psf_fwhm_area'], marker='.', s=1, alpha=0.01)\n", + "plt.hexbin()\n", + "plt.plot(np.linspace(0.5, 1.5, 10), np.linspace(0.5, 1.5, 10), c='black', ls='--')\n", + "plt.plot(np.linspace(0.5, 1.5, 10), np.linspace(0.5, 1.5, 10) + 0.12, c='red', ls='--')\n", + "plt.xlim(0.5, 1.5)\n", + "plt.ylim(0.5, 1.5)\n", + "plt.xlabel('PSF (moments)')\n", + "plt.ylabel('PSF (area)')\n", + "plt.grid()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "68c370cb-989a-4be0-a2b6-bfaa36638d19", + "metadata": {}, + "outputs": [], + "source": [ + "ccdvisits['target_name'].str.contains('ddf')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "268c1f5b-eb2a-41a8-b1c2-9ae84401dea1", + "metadata": {}, + "outputs": [], + "source": [ + "bins = np.linspace(0.5, 3., 50)\n", + "\n", + "selection = ccdvisits['target_name'].str.contains('ddf') & (np.fabs(ccdvisits['physical_rotator_angle']) < 20.)\n", + "selection = ccdvisits['target_name'].str.contains('ddf') & (np.fabs(ccdvisits['physical_rotator_angle']) > 20.)\n", + "#selection_2 = ccdvisits['target_name'].str.contains('ddf') & (ccdvisits['airmass'] > 1.2)\n", + "\n", + "#print(np.sum(selection) / 181)\n", + "#print(np.sum(selection_2) / 181)\n", + "\n", + "print(np.unique(ccdvisits['target_name'][selection]))\n", + "\n", + "plt.figure()\n", + "plt.hist(ccdvisits[\"psf_fwhm\"][selection], bins=bins, histtype='step', density=True, label='DDF')\n", + "plt.hist(ccdvisits[\"psf_fwhm\"][selection_2], bins=bins, histtype='step', density=True, label='DDF')\n", + "#plt.hist(ccdvisits[\"psf_fwhm\"][~selection], bins=bins, histtype='step', density=True, label='non-DDF')\n", + "plt.legend()\n", + "plt.title('FBS Visits: %s - %s'%(day_obs_min, day_obs_max))\n", + "plt.xlabel('PSF FWHM (arcsec)')\n", + "plt.ylabel('PDF')\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "70226b1d-eb2b-42da-87da-bdce5c2595af", + "metadata": {}, + "outputs": [], + "source": [ + "selection = ccdvisits['target_name'].str.contains('ddf')\n", + "\n", + "plt.figure()\n", + "for band in bands:\n", + " band_selection = selection & (ccdvisits[\"band\"] == band)\n", + " plt.scatter(\n", + " ccdvisits[\"seq_num\"][selection & band_selection], \n", + " ccdvisits[\"psf_fwhm\"][selection & band_selection],\n", + " c=colors[band], label=band,\n", + " edgecolor='none', s=1\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "63e301a3-34f9-4d3a-9fe5-f6efee445e44", + "metadata": {}, + "outputs": [], + "source": [ + "ccdvisits" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c34f7960-f1e9-4e12-bb96-8b770a303230", + "metadata": {}, + "outputs": [], + "source": [ + "np.sort(ccdvisits['seq_num'].unique())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "087e45dc-af16-46b3-aee9-5f18806eefa7", + "metadata": {}, + "outputs": [], + "source": [ + "selection = (ccdvisits['band'] != \"y\") & (ccdvisits['day_obs'] >= 20251026)\n", + "print(np.sum(selection))\n", + "\n", + "plt.figure()\n", + "plt.hexbin(\n", + " ccdvisits[\"donut_blur_atm_fwhm\"][selection], \n", + " ccdvisits[\"psf_fwhm\"][selection], \n", + " #ccdvisits[\"psf_fwhm_zenith_500nm\"][selection], \n", + " gridsize=50, mincnt=10, extent=(0.5, 2.0, 0.5, 2.0), cmap='viridis')\n", + "plt.colorbar(label='Density')\n", + "plt.plot(np.linspace(0.5, 2.0, 10), np.linspace(0.5, 2.0, 10), lw=2, ls='--', c='black')\n", + "\n", + "x = np.linspace(0.5, 2.0, 100)\n", + "plt.plot(x, np.sqrt(x**2 + 0.4**2), c='red', lw=2, ls='--', label='0.4\" System Contribution')\n", + "plt.plot(x, np.sqrt(x**2 + 0.5**2), c='orange', lw=2, ls='--', label='0.5\" System Contribution')\n", + "plt.plot(x, np.sqrt(x**2 + 0.6**2), c='cyan', lw=2, ls='--', label='0.6\" System Contribution')\n", + "plt.plot(x, np.sqrt(x**2 + 0.7**2), c='magenta', lw=2, ls='--', label='0.7\" System Contribution')\n", + "\n", + "plt.xlim(0.5, 1.8)\n", + "plt.ylim(0.5, 1.8)\n", + "\n", + "plt.legend(loc='lower right')\n", + "\n", + "#plt.xlabel(\"$\\sqrt{ {\\rm Donut_Blur_FWHM}^2 - {\\rm CAM_FWHM}^2}$ (arcsec)\")\n", + "plt.xlabel(\"$\\sqrt{\\mathrm{DONUT\\_BLUR\\_FWHM}^2 - \\mathrm{CAM\\_FWHM}^2}$ (arcsec)\")\n", + "#plt.ylabel(\"psf_fwhm_zenith_500nm\")\n", + "plt.ylabel(\"PSF FWHM (arcsec)\")\n", + "plt.title(\"Using Donut Blur as Proxy for Atmosphere Contribution\")\n", + "#plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "281a2ef5-cf1d-43a7-ba77-8219bad631f5", + "metadata": {}, + "outputs": [], + "source": [ + "bins = np.linspace(0., 2., 41)\n", + "plt.figure()\n", + "plt.hist(\n", + " np.sqrt(ccdvisits[\"psf_fwhm\"][selection]**2 + CAM_FWHM**2 - ccdvisits[\"donut_blur_fwhm\"][selection]**2), \n", + " bins=bins, histtype='step', label='w/ CAM_FWHM',\n", + ")\n", + "plt.hist(\n", + " np.sqrt(ccdvisits[\"psf_fwhm\"][selection]**2 - ccdvisits[\"donut_blur_fwhm\"][selection]**2), \n", + " bins=bins, histtype='step', label='w/out CAM_FWHM',\n", + ")\n", + "plt.legend()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a4c050f2-2873-4c3b-879e-371a2bf7ce04", + "metadata": {}, + "outputs": [], + "source": [ + "bins = np.linspace(0., 2., 41)\n", + "plt.figure()\n", + "plt.hist(\n", + " ccdvisits[\"aos_fwhm\"][selection], \n", + " bins=bins, histtype='step', label='w/out CAM_FWHM',\n", + ")\n", + "plt.hist(\n", + " ccdvisits[\"aos_cam_fwhm\"], \n", + " bins=bins, histtype='step', label='w/ CAM_FWHM',\n", + ")\n", + "plt.legend()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b04a8492-1772-473d-a32d-63a6bcd9fc26", + "metadata": {}, + "outputs": [], + "source": [ + "\"\"\"\n", + "selection = (ccdvisits['seq_num'] == 52) \n", + "\n", + "bins = np.linspace(-0.5, 0.5, 150)\n", + "plt.figure()\n", + "plt.hist(ccdvisits['psf_e1'][selection], bins=bins, histtype='step', label='e1')\n", + "plt.hist(ccdvisits['psf_e2'][selection], bins=bins, histtype='step', label='e2')\n", + "plt.hist(ccdvisits['psf_e'][selection], bins=bins, histtype='step', label='e')\n", + "plt.legend()\n", + "\"\"\"" + ] + }, + { + "cell_type": "markdown", + "id": "a63c7e92-7c44-415e-bce9-1ca30fe90703", + "metadata": {}, + "source": [ + "## Publication Style Plots" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "60e8436d-36e5-403b-a6b3-1d343457f6ad", + "metadata": {}, + "outputs": [], + "source": [ + "publication_plots.set_rubin_plotstyle()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f1928fa1-85a6-420c-8dfa-fd274c044d57", + "metadata": {}, + "outputs": [], + "source": [ + "#selection = (ccdvisits['day_obs'] >= 20251130) & (ccdvisits['day_obs'] <= 20251202)\n", + "selection = (ccdvisits['day_obs'] >= day_obs_min) & (ccdvisits['day_obs'] <= day_obs_max)\n", + "print(np.sum(selection) / 189.)\n", + "\n", + "\n", + "bins = np.linspace(-0.4, 0.4, 81)\n", + "\n", + "plt.figure()\n", + "plt.hist(ccdvisits['psf_e'][selection], bins=bins, density=True, histtype='step', lw=2, label='PSF e')\n", + "#plt.hist(ccdvisits['psf_e1'][selection], bins=bins, density=True, histtype='step', lw=2, label='PSF e1')\n", + "#plt.hist(ccdvisits['psf_e2'][selection], bins=bins, density=True, histtype='step', lw=2, label='PSF e2')\n", + "plt.legend()\n", + "plt.xlabel('Ellipticity')\n", + "plt.ylabel('PDF')\n", + "#plt.xlim(bins[0], bins[-1])\n", + "plt.xlim(0, bins[-1])\n", + "plt.ylim(0., plt.ylim()[-1])\n", + "plt.title('%i - %s'%(min(ccdvisits['day_obs'][selection]), max(ccdvisits['day_obs'][selection])))\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "961a79e9-84ae-40e6-8705-2e7db94ee0d3", + "metadata": {}, + "outputs": [], + "source": [ + "bins = np.linspace(0.5, 2., 51)\n", + "\n", + "plt.figure()\n", + "\n", + "for band in bands:\n", + " selection = (ccdvisits[\"band\"] == band)\n", + " plt.hist(ccdvisits[\"psf_fwhm\"][selection], bins=bins, color=colors[band], density=True, histtype=\"step\", lw=2, label=band)\n", + "\n", + "plt.grid()\n", + "plt.xlabel('PSF FWHM (arcsec)')\n", + "plt.ylabel('Fraction of Sensors')\n", + "plt.legend()\n", + "plt.xlim(np.min(bins), np.max(bins))\n", + "plt.title(f'{day_obs_min}-{day_obs_max}')\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ba7647ba-22a4-4634-9ded-40d256e000fd", + "metadata": {}, + "outputs": [], + "source": [ + "bins = np.linspace(0.5, 2.5, 51)\n", + "\n", + "plt.figure()\n", + "\n", + "for band in bands:\n", + " selection = (ccdvisits[\"band\"] == band)\n", + " plt.hist(ccdvisits[\"psf_fwhm\"][selection], bins=bins, color=colors[band], density=True, histtype=\"step\", cumulative=True, lw=2, label=band)\n", + "plt.grid()\n", + "plt.xlabel('PSF FWHM (arcsec)')\n", + "plt.ylabel('Fraction of Sensors')\n", + "plt.legend(loc='lower right')\n", + "plt.xlim(np.min(bins), np.max(bins))\n", + "plt.ylim(0., 1.)\n", + "plt.title(f'{day_obs_min}-{day_obs_max}')\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f71b2746-97c1-422e-9133-f1507b396c17", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "groups = ccdvisits.groupby('visit_id')\n", + "visits_summary = pd.DataFrame({\n", + " 'day_obs': groups['day_obs'].first(),\n", + " 'target_name': groups['target_name'].first(),\n", + " 'science_program': groups['science_program'].first(),\n", + " 'observation_reason': groups['observation_reason'].first(),\n", + " 'seq_num': groups['seq_num'].median(),\n", + " 'exp_midpt_mjd': groups['exp_midpt_mjd'].median(),\n", + " 'donut_blur_fwhm': groups['donut_blur_fwhm'].median(),\n", + " 'aos_fwhm': groups['aos_fwhm'].median(),\n", + " 'donut_blur_atm_fwhm': groups['donut_blur_atm_fwhm'].median(),\n", + " 'aos_cam_fwhm': groups['aos_cam_fwhm'].median(),\n", + " 'physical_rotator_angle': groups['physical_rotator_angle'].median(),\n", + " 'altitude': groups['altitude'].median(),\n", + " 'airmass': groups['airmass'].median(),\n", + " 'psf_fwhm_05': groups['psf_fwhm'].quantile(0.05),\n", + " 'psf_fwhm_50': groups['psf_fwhm'].quantile(0.50),\n", + " 'psf_fwhm_95': groups['psf_fwhm'].quantile(0.95),\n", + " 'psf_fwhm_zenith_500nm_50': groups['psf_fwhm_zenith_500nm'].quantile(0.50),\n", + " 'psf_fwhm_area_50': groups['psf_fwhm_area'].quantile(0.50),\n", + " 'psf_e_05': groups['psf_e'].quantile(0.05),\n", + " 'psf_e_50': groups['psf_e'].quantile(0.50),\n", + " 'psf_e_95': groups['psf_e'].quantile(0.95),\n", + " 'psf_e1_05': groups['psf_e1'].quantile(0.05),\n", + " 'psf_e1_50': groups['psf_e1'].quantile(0.50),\n", + " 'psf_e1_95': groups['psf_e1'].quantile(0.95),\n", + " 'psf_e2_05': groups['psf_e2'].quantile(0.05),\n", + " 'psf_e2_50': groups['psf_e2'].quantile(0.50),\n", + " 'psf_e2_95': groups['psf_e2'].quantile(0.95),\n", + " 'psf_e_unnormalized_05': groups['psf_e_unnormalized'].quantile(0.05),\n", + " 'psf_e_unnormalized_50': groups['psf_e_unnormalized'].quantile(0.50),\n", + " 'psf_e_unnormalized_95': groups['psf_e_unnormalized'].quantile(0.95),\n", + " 'band': groups['band'].first(),\n", + "})\n", + "visits_summary['psf_fwhm_95_05'] = np.sqrt(visits_summary['psf_fwhm_95']**2 - visits_summary['psf_fwhm_05']**2)\n", + "visits_summary['sys_50'] = np.sqrt(visits_summary['psf_fwhm_50']**2 + CAM_FWHM**2 - visits_summary['donut_blur_fwhm']**2)\n", + "visits_summary['psf_fwhm_model'] = np.sqrt(visits_summary['aos_fwhm']**2 + visits_summary['donut_blur_fwhm']**2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1eb8f692-d6ec-433f-924b-430529ef6c65", + "metadata": {}, + "outputs": [], + "source": [ + "bins = np.linspace(0.2, 1.0, 51)\n", + "\n", + "plt.figure()\n", + "\n", + "plt.hist(visits_summary[\"aos_cam_fwhm\"], bins=bins, color='tab:blue', density=True, histtype=\"step\", cumulative=True, lw=2)\n", + "plt.grid()\n", + "plt.xlabel('Instrument Contribution to Image Quality:\\nOptics + Camera Diffusion (arcsec)')\n", + "plt.ylabel('Fraction of Visits')\n", + "plt.xlim(np.min(bins), np.max(bins))\n", + "plt.ylim(0., 1.)\n", + "plt.title(f'{day_obs_min}-{day_obs_max}')\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cf9ebd1e-3e2d-4bd0-a79d-e16b3201cda4", + "metadata": {}, + "outputs": [], + "source": [ + "visits_summary['observation_reason']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d059d6be-aa54-49f6-830b-be2ef28270c7", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "plt.scatter(visits_summary['psf_fwhm_50'], visits_summary['psf_e_95'], s=2)\n", + "#plt.scatter(visits_summary['donut_blur_fwhm'], visits_summary['psf_e_50'], s=2)\n", + "#plt.scatter(visits_summary['psf_fwhm_95_05'], visits_summary['psf_e_50'], s=2)\n", + "plt.axvline(0.8, c='0.75', ls='--')\n", + "plt.axhline(0.07, c='0.75', ls='--')\n", + "plt.xlabel('Per-visit Median PSF FWHM (arcsec)')\n", + "#plt.ylabel('Per-visit Mediain PSF e')\n", + "plt.ylabel('Per-visit 95th Percentile PSF e')\n", + "plt.title('%i - %i'%(day_obs_min, day_obs_max))\n", + "plt.ylim(0., 0.3)\n", + "plt.xlim(0.65, 1.2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6a40491a-0f51-4a42-bda6-21ef991b14ab", + "metadata": {}, + "outputs": [], + "source": [ + "np.median(visits_summary['psf_e_95'])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6a12688c-6292-4d5f-aa4b-cf0007feaae4", + "metadata": {}, + "outputs": [], + "source": [ + "visits_summary" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0e97f8e0-b897-4ded-93a5-b54136ecd1e0", + "metadata": {}, + "outputs": [], + "source": [ + "print(visits_summary['day_obs'].value_counts().sort_index().to_string())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ae80048d-e8b9-4405-9c5a-2cad0a0cce1e", + "metadata": {}, + "outputs": [], + "source": [ + "np.mean(visits_summary['day_obs'].value_counts()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bc99fc6e-5d00-4ef1-b5fd-6aed5a951ac0", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "plt.scatter(visits_summary['aos_fwhm'], visits_summary['psf_fwhm_95_05'], s=1)\n", + "#plt.xlim(0., 1.)\n", + "#plt.ylim(0., 1.)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4260aaf1-adf5-42fa-8299-26d89f79f50e", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "plt.hist(visits_summary['psf_fwhm_50'], bins=41)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "427d4c1b-063f-4e06-be52-c127326c68bb", + "metadata": {}, + "outputs": [], + "source": [ + "selection = (visits_summary['band'] != \"y\") & (visits_summary['day_obs'] >= 20251026)\n", + "print(np.sum(selection))\n", + "\n", + "plt.figure()\n", + "plt.hexbin(\n", + " visits_summary[\"donut_blur_atm_fwhm\"][selection], \n", + " visits_summary[\"psf_fwhm_50\"][selection], \n", + " #visits_summary[\"psf_fwhm_zenith_500nm_50\"][selection], \n", + " gridsize=40, mincnt=1, extent=(0.5, 2.0, 0.5, 2.0), cmap='viridis')\n", + "plt.colorbar(label='Density')\n", + "plt.plot(np.linspace(0.5, 2.0, 10), np.linspace(0.5, 2.0, 10), lw=2, ls='--', c='black')\n", + "\n", + "x = np.linspace(0.5, 2.0, 100)\n", + "plt.plot(x, np.sqrt(x**2 + 0.4**2), c='red', lw=2, ls='--', label='0.4\" System Contribution')\n", + "plt.plot(x, np.sqrt(x**2 + 0.5**2), c='orange', lw=2, ls='--', label='0.5\" System Contribution')\n", + "plt.plot(x, np.sqrt(x**2 + 0.6**2), c='cyan', lw=2, ls='--', label='0.6\" System Contribution')\n", + "plt.plot(x, np.sqrt(x**2 + 0.7**2), c='magenta', lw=2, ls='--', label='0.7\" System Contribution')\n", + "\n", + "plt.xlim(0.5, 1.8)\n", + "plt.ylim(0.5, 1.8)\n", + "\n", + "#plt.xlabel(\"Donut Blur FWHM (arcsec)\")\n", + "plt.xlabel(\"$\\sqrt{\\mathrm{DONUT\\_BLUR\\_FWHM}^2 - \\mathrm{CAM\\_FWHM}^2}$ (arcsec)\")\n", + "#plt.ylabel(\"psf_fwhm_zenith_500nm\")\n", + "plt.ylabel(\"PSF FWHM (arcsec)\")\n", + "plt.title(\"Using Donut Blur as Proxy for Atmosphere Contribution\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "b63ddde4-3e53-4f1d-b9b0-8aefaedae372", + "metadata": { + "jp-MarkdownHeadingCollapsed": true + }, + "source": [ + "## Comparison to Simulation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "61a9411d-f460-4e99-9ef5-c2402e415e64", + "metadata": {}, + "outputs": [], + "source": [ + "import sqlite3" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e7ac7b67-a907-4260-b7af-5b461f79262f", + "metadata": {}, + "outputs": [], + "source": [ + "db = '/sdf/group/rubin/web_data/sim-data/sims_featureScheduler_runs5.1/baseline/baseline_v5.1.1_10yrs.db'\n", + "\n", + "conn = sqlite3.connect(db)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ca30b5cd-171f-407e-9623-47974160258d", + "metadata": {}, + "outputs": [], + "source": [ + "observations = pd.read_sql('select * from observations;', conn)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5f991a88-867a-47f8-9534-3a2689cdebcc", + "metadata": {}, + "outputs": [], + "source": [ + "selection_season = (observations['night'] == 0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7c494ef6-68c6-4c2a-a36a-0176ddce4039", + "metadata": {}, + "outputs": [], + "source": [ + "from astropy.time import Time\n", + "\n", + "time_start_sim = Time(observations['observationStartMJD'][selection_season][0], format='mjd')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dce723d9-91c9-478b-b569-a1d8f73c2fb6", + "metadata": {}, + "outputs": [], + "source": [ + "def dayObsToTime(day_obs):\n", + " year = str(day_obs)[0:4]\n", + " month = str(day_obs)[4:6]\n", + " day = str(day_obs)[6:8]\n", + " time = Time('%s-%s-%s 00:00:00'%(year, month, day), format='iso')\n", + " return time" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5d07670a-b745-47b5-8573-5d490e06cf2e", + "metadata": {}, + "outputs": [], + "source": [ + "time_day_obs_min = dayObsToTime(day_obs_min)\n", + "time_day_obs_max = dayObsToTime(day_obs_max)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dd8a43bc-7833-428c-aa46-4361320cd418", + "metadata": {}, + "outputs": [], + "source": [ + "(time_day_obs_min - time_start_sim).value" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "51585d1d-91df-4eb5-a5c9-c56cdf971432", + "metadata": {}, + "outputs": [], + "source": [ + "selection_season = ((observations['night'] % 365) > (time_day_obs_min - time_start_sim).value) & ((observations['night'] % 365) < (time_day_obs_max - time_start_sim).value) " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4ae782db-06b9-4810-96cd-66c8106b4ae2", + "metadata": {}, + "outputs": [], + "source": [ + "Time(observations['observationStartMJD'][selection_season], format='mjd')[-1].iso" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ce7bb7e8-f4ab-4fc7-8142-c0ce875847c0", + "metadata": {}, + "outputs": [], + "source": [ + "bins = np.linspace(0.5, 2., 401)\n", + "\n", + "plt.figure()\n", + "\n", + "for band in bands:\n", + " selection = (observations[\"band\"] == band)\n", + " plt.hist(observations['seeingFwhmEff'][selection], bins=bins, color=colors[band], density=True, histtype=\"step\", cumulative=True, label=band, alpha=0.5)\n", + "\n", + " selection = (ccdvisits[\"band\"] == band)\n", + " plt.hist(ccdvisits[\"psf_fwhm_area\"][selection], bins=bins, color=colors[band], density=True, histtype=\"step\", cumulative=True, label=band, lw=2)\n", + "plt.grid()\n", + "plt.xlabel('PSF FWHM from area (arcsec)')\n", + "plt.ylabel('Fraction of Visits')\n", + "plt.legend(loc='lower right')\n", + "plt.xlim(np.min(bins), np.max(bins))\n", + "plt.ylim(0., 1.)\n", + "#plt.title(f'{day_obs_min}-{day_obs_max}')\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "67afc00b-78f9-4f1d-aa05-c26b892520e3", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(9, 6), nrows=2, ncols=3)\n", + "\n", + "bins = np.linspace(0.5, 2., 401)\n", + "\n", + "airmass_min = 1.2\n", + "airmass_max = 1.4\n", + "\n", + "for ii, band in enumerate(bands):\n", + " row, column = int(ii/3), ii%3\n", + " \n", + " selection = (observations[\"band\"] == band)\n", + " selection_airmass = (observations[\"airmass\"] > airmass_min) & (observations[\"airmass\"] < airmass_max)\n", + " #quantity = observations['seeingFwhmEff'][selection & selection_season & selection_airmass]\n", + " quantity = observations['seeingFwhmGeom'][selection & selection_season & selection_airmass]\n", + " ax[row][column].hist(quantity, bins=bins, color=colors[band], density=True, histtype=\"step\", cumulative=True, \n", + " label=f'{band} sim\\nmedian=%.2f'%(np.median(quantity)), alpha=0.5)\n", + "\n", + " #selection = (ccdvisits[\"band\"] == band)\n", + " #ax[row][column].hist(ccdvisits[\"psf_fwhm\"][selection], bins=bins, color=colors[band], density=True, histtype=\"step\", cumulative=True, label=f'{band} measured', lw=2)\n", + " selection = (visits_summary[\"band\"] == band)\n", + " selection_airmass = (visits_summary[\"airmass\"] > airmass_min) & (visits_summary[\"airmass\"] < airmass_max)\n", + " #ax[row][column].hist(visits_summary[\"psf_fwhm_50\"][selection], bins=bins, color=colors[band], density=True, histtype=\"step\", cumulative=True, label=f'{band} measured (moments)', lw=2, alpha=0.75)\n", + " quantity = visits_summary[\"psf_fwhm_50\"][selection & selection_airmass]\n", + " #quantity = visits_summary[\"psf_fwhm_area_50\"][selection & selection_airmass]\n", + " ax[row][column].hist(quantity, bins=bins, color=colors[band], density=True, histtype=\"step\", cumulative=True, \n", + " label=f'{band} measured\\nmedian=%.2f'%(np.median(quantity)), lw=2)\n", + "\n", + " ax[row][column].grid()\n", + " ax[row][column].set_xlabel('PSF FWHM (arcsec)')\n", + " ax[row][column].set_ylabel('Fraction of Visits')\n", + " ax[row][column].legend(loc='lower right', fontsize=10)\n", + " ax[row][column].set_xlim(np.min(bins), np.max(bins))\n", + " ax[row][column].set_ylim(0., 1.)\n", + "#plt.title(f'{day_obs_min}-{day_obs_max}')\n", + "plt.suptitle(f'{day_obs_min} - {day_obs_max}; {airmass_min} < airmass < {airmass_max}')\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f14e225f-c627-4546-87ce-1a2a6044b713", + "metadata": {}, + "outputs": [], + "source": [ + "for ii in range(0, 6):\n", + " print(int(ii/3), ii%3)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7d7561e6-4cec-43ec-85b5-3946a6606e88", + "metadata": {}, + "outputs": [], + "source": [ + "np.sqrt(1.0**2 + 0.7**2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "45391b12-fa65-4d7a-a44e-3455295dbd1d", + "metadata": {}, + "outputs": [], + "source": [ + "#import matplotlib.cm as cm\n", + "#cmap = cm.get_cmap('viridis')\n", + "\n", + "plt.figure(figsize=(6,6))\n", + "\n", + "#instrument = 0.5\n", + "\n", + "for tot in [0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2]:\n", + " sys = np.linspace(0.3, 1.5, 1000)\n", + " atm = np.sqrt(tot**2 - sys**2)\n", + " plt.plot(sys, atm, label='%.1f'%(tot), lw=3)\n", + " #dome = np.linspace(0.2, 1.5, 1000)\n", + " #atm = np.sqrt(tot**2 - (instrument**2 + dome**2))\n", + " #plt.plot(dome, atm, label='%.1f'%(tot))\n", + "\n", + "plt.xlabel('System Contribution (arscec)')\n", + "plt.ylabel('Atmosphere Contribution (arcsec)')\n", + "plt.xlim(0.3, 1.2)\n", + "plt.ylim(0.3, 1.2)\n", + "plt.legend(title='Total DIQ (arcsec)')\n", + "plt.grid()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f069502c-d906-4b5a-a1ed-1961ac1d1085", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(figsize=(6,6))\n", + "\n", + "for atm in [0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1]:\n", + " sys = np.linspace(0.3, 1.5, 1000)\n", + " tot = np.sqrt(atm**2 + sys**2)\n", + " plt.plot(sys, tot, label='%.1f'%(atm), lw=3)\n", + " #dome = np.linspace(0.2, 1.5, 1000)\n", + " #atm = np.sqrt(tot**2 - (instrument**2 + dome**2))\n", + " #plt.plot(dome, atm, label='%.1f'%(tot))\n", + "\n", + "plt.xlabel('System Contribution (arscec)')\n", + "plt.ylabel('Total DIQ (arcsec)')\n", + "plt.xlim(0.3, 1.2)\n", + "plt.ylim(0.5, 1.2)\n", + "plt.legend(title='Atmosphere Contribution (arcsec)')\n", + "plt.grid()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f2f2a853-1f61-458a-9976-b91335bef420", + "metadata": {}, + "outputs": [], + "source": [ + "np.sqrt(0.5**2 + 0.2**2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "00b2f458-f45e-4d70-a406-0c10b60876e6", + "metadata": {}, + "outputs": [], + "source": [ + "bins = np.linspace(1., 2.2, 41)\n", + "\n", + "plt.figure()\n", + "plt.hist(visits_summary[\"airmass\"], bins=bins, density=True, histtype='step', label='measured')\n", + "plt.hist(observations['airmass'][selection_season], bins=bins, density=True, histtype='step', label='sim')\n", + "plt.legend()\n", + "plt.xlabel('Airmass')\n", + "plt.ylabel('PDF')\n", + "plt.title(f'{day_obs_min} - {day_obs_max}')" + ] + }, + { + "cell_type": "markdown", + "id": "0a8c8f9e-5177-4bad-b80a-46495cec9e7d", + "metadata": { + "jp-MarkdownHeadingCollapsed": true + }, + "source": [ + "## Zernikes" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a438df67-83c3-4ef6-a52e-938083a44738", + "metadata": {}, + "outputs": [], + "source": [ + "import galsim\n", + "\n", + "def getPsfGradPerZernike(\n", + " diameter: float = 8.36,\n", + " obscuration: float = 0.612,\n", + " jmin: int = 4,\n", + " jmax: int = 22,\n", + ") -> np.ndarray:\n", + " \"\"\"Get the gradient of the PSF FWHM with respect to each Zernike.\n", + "\n", + " This function takes no positional arguments. All parameters must be passed\n", + " by name (see the list of parameters below).\n", + "\n", + " Parameters\n", + " ----------\n", + " diameter : float, optional\n", + " The diameter of the telescope aperture, in meters.\n", + " (the default, 8.36, corresponds to the LSST primary mirror)\n", + " obscuration : float, optional\n", + " Central obscuration of telescope aperture (i.e. R_outer / R_inner).\n", + " (the default, 0.612, corresponds to the LSST primary mirror)\n", + " jmin : int, optional\n", + " The minimum Noll index, inclusive. Must be >= 0. (the default is 4)\n", + " jmax : int, optional\n", + " The max Zernike Noll index, inclusive. Must be >= jmin.\n", + " (the default is 22.)\n", + "\n", + " Returns\n", + " -------\n", + " np.ndarray\n", + " Gradient of the PSF FWHM with respect to the corresponding Zernike.\n", + " Units are arcsec / micron.\n", + "\n", + " Raises\n", + " ------\n", + " ValueError\n", + " If jmin is negative or jmax is less than jmin\n", + " \"\"\"\n", + " # Check jmin and jmax\n", + " if jmin < 0:\n", + " raise ValueError(\"jmin cannot be negative.\")\n", + " if jmax < jmin:\n", + " raise ValueError(\"jmax must be greater than jmin.\")\n", + "\n", + " # Calculate the conversion factors\n", + " conversion_factors = np.zeros(jmax + 1)\n", + " for i in range(jmin, jmax + 1):\n", + " # Set coefficients for this Noll index: coefs = [0, 0, ..., 1]\n", + " # Note the first coefficient is Noll index 0, which does not exist and\n", + " # is therefore always ignored by galsim\n", + " coefs = [0] * i + [1]\n", + "\n", + " # Create the Zernike polynomial with these coefficients\n", + " R_outer = diameter / 2\n", + " R_inner = R_outer * obscuration\n", + " Z = galsim.zernike.Zernike(coefs, R_outer=R_outer, R_inner=R_inner)\n", + "\n", + " # We can calculate the size of the PSF from the RMS of the gradient of\n", + " # the wavefront. The gradient of the wavefront perturbs photon paths.\n", + " # The RMS quantifies the size of the collective perturbation.\n", + " # If we expand the wavefront gradient in another series of Zernike\n", + " # polynomials, we can exploit the orthonormality of the Zernikes to\n", + " # calculate the RMS from the Zernike coefficients.\n", + " rms_tilt = np.sqrt(np.sum(Z.gradX.coef**2 + Z.gradY.coef**2) / 2)\n", + "\n", + " # Convert to arcsec per micron\n", + " rms_tilt = np.rad2deg(rms_tilt * 1e-6) * 3600\n", + "\n", + " # Convert rms -> fwhm\n", + " fwhm_tilt = 2 * np.sqrt(2 * np.log(2)) * rms_tilt\n", + "\n", + " # Save this conversion factor\n", + " conversion_factors[i] = fwhm_tilt\n", + "\n", + " return conversion_factors[jmin:]\n", + "\n", + "\n", + "def convertZernikesToPsfWidth(\n", + " zernikes: np.ndarray,\n", + " diameter: float = 8.36,\n", + " obscuration: float = 0.612,\n", + " jmin: int = 4,\n", + ") -> np.ndarray:\n", + " \"\"\"Convert Zernike amplitudes to quadrature contribution to the PSF FWHM.\n", + "\n", + " Parameters\n", + " ----------\n", + " zernikes : np.ndarray\n", + " Zernike amplitudes (in microns), starting with Noll index `jmin`.\n", + " Either a 1D array of zernike amplitudes, or a 2D array, where each row\n", + " corresponds to a different set of amplitudes.\n", + " diameter : float\n", + " The diameter of the telescope aperture, in meters.\n", + " (the default, 8.36, corresponds to the LSST primary mirror)\n", + " obscuration : float\n", + " Central obscuration of telescope aperture (i.e. R_outer / R_inner).\n", + " (the default, 0.612, corresponds to the LSST primary mirror)\n", + " jmin : int\n", + " The minimum Zernike Noll index, inclusive. Must be >= 0. The\n", + " max Noll index is inferred from `jmin` and the length of `zernikes`.\n", + " (the default is 4, which ignores piston, x & y offsets, and tilt.)\n", + "\n", + " Returns\n", + " -------\n", + " dFWHM: np.ndarray\n", + " Quadrature contribution of each Zernike vector to the PSF FWHM\n", + " (in arcseconds).\n", + "\n", + " Notes\n", + " -----\n", + " Converting Zernike amplitudes to their quadrature contributions to the PSF\n", + " FWHM allows for easier physical interpretation of Zernike amplitudes and\n", + " the performance of the AOS system.\n", + "\n", + " For example, image we have a true set of zernikes, [Z4, Z5, Z6], such that\n", + " ConvertZernikesToPsfWidth([Z4, Z5, Z6]) = [0.1, -0.2, 0.3] arcsecs.\n", + " These Zernike perturbations increase the PSF FWHM by\n", + " sqrt[(0.1)^2 + (-0.2)^2 + (0.3)^2] ~ 0.37 arcsecs.\n", + "\n", + " If the AOS perfectly corrects for these perturbations, the PSF FWHM will\n", + " not increase in size. However, imagine the AOS estimates zernikes, such\n", + " that ConvertZernikesToPsfWidth([Z4, Z5, Z6]) = [0.1, -0.3, 0.4] arcsecs.\n", + " These estimated Zernikes, do not exactly match the true Zernikes above.\n", + " Therefore, the post-correction PSF will still be degraded with respect to\n", + " the optimal PSF. In particular, the PSF FWHM will be increased by\n", + " sqrt[(0.1 - 0.1)^2 + (-0.2 - (-0.3))^2 + (0.3 - 0.4)^2] ~ 0.14 arcsecs.\n", + "\n", + " This conversion depends on a linear approximation that begins to break down\n", + " for RSS(dFWHM) > 0.20 arcsecs. Beyond this point, the approximation tends\n", + " to overestimate the PSF degradation. In other words, if\n", + " sqrt(sum( dFWHM^2 )) > 0.20 arcsec, it is likely that dFWHM is\n", + " over-estimated. However, the point beyond which this breakdown begins\n", + " (and whether the approximation over- or under-estimates dFWHM) can change,\n", + " depending on which Zernikes have large amplitudes. In general, if you have\n", + " large Zernike amplitudes, proceed with caution!\n", + " Note that if the amplitudes Z_est and Z_true are large, this is okay, as\n", + " long as |Z_est - Z_true| is small.\n", + "\n", + " For a notebook demonstrating where the approximation breaks down:\n", + " https://gist.github.com/jfcrenshaw/24056516cfa3ce0237e39507674a43e1\n", + "\n", + " Raises\n", + " ------\n", + " ValueError\n", + " If jmin is negative\n", + " \"\"\"\n", + " # Check jmin\n", + " if jmin < 0:\n", + " raise ValueError(\"jmin cannot be negative.\")\n", + "\n", + " # Calculate jmax from jmin and the length of the zernike array\n", + " jmax = jmin + np.array(zernikes).shape[-1] - 1\n", + "\n", + " # Calculate the conversion factors for each zernike\n", + " conversion_factors = getPsfGradPerZernike(\n", + " jmin=jmin,\n", + " jmax=jmax,\n", + " diameter=diameter,\n", + " obscuration=obscuration,\n", + " )\n", + "\n", + " # Convert the Zernike amplitudes from microns to their quadrature\n", + " # contribution to the PSF FWHM\n", + " dFWHM = conversion_factors * zernikes\n", + "\n", + " return dFWHM" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fa2b3ba7-6328-419e-880b-0ccff339d3c7", + "metadata": {}, + "outputs": [], + "source": [ + "def fetchZernikes(day_obs_min, day_obs_max):\n", + " \n", + " visits_query = f'''\n", + " SELECT \n", + " ccdvisit1_quicklook.z4,\n", + " ccdvisit1_quicklook.z5,\n", + " ccdvisit1_quicklook.z6,\n", + " ccdvisit1_quicklook.z7,\n", + " ccdvisit1_quicklook.z8,\n", + " ccdvisit1_quicklook.z9,\n", + " ccdvisit1_quicklook.z10,\n", + " ccdvisit1_quicklook.z11,\n", + " ccdvisit1_quicklook.z12,\n", + " ccdvisit1_quicklook.z13,\n", + " ccdvisit1_quicklook.z14,\n", + " ccdvisit1_quicklook.z15,\n", + " ccdvisit1_quicklook.z16,\n", + " ccdvisit1_quicklook.z17,\n", + " ccdvisit1_quicklook.z18,\n", + " ccdvisit1_quicklook.z19,\n", + " ccdvisit1_quicklook.z20,\n", + " ccdvisit1_quicklook.z21,\n", + " ccdvisit1_quicklook.z22,\n", + " ccdvisit1_quicklook.z23,\n", + " ccdvisit1_quicklook.z24,\n", + " ccdvisit1_quicklook.z25,\n", + " ccdvisit1_quicklook.z26,\n", + " ccdvisit1_quicklook.z27,\n", + " ccdvisit1_quicklook.z28,\n", + " ccdvisit1.detector,\n", + " visit1.visit_id\n", + " FROM\n", + " cdb_{instrument}.ccdvisit1_quicklook AS ccdvisit1_quicklook,\n", + " cdb_{instrument}.ccdvisit1 AS ccdvisit1,\n", + " cdb_{instrument}.visit1 AS visit1 \n", + " WHERE \n", + " ccdvisit1.ccdvisit_id = ccdvisit1_quicklook.ccdvisit_id\n", + " AND ccdvisit1.visit_id = visit1.visit_id \n", + " AND ccdvisit1.detector IN (191, 195, 199, 200, 203) -- corner wavefront sensors\n", + " -- AND ccdvisit1.detector IN (191, 192, 195, 196, 199, 200, 203, 204) -- corner wavefront sensors\n", + " AND visit1.airmass > 0\n", + " AND visit1.day_obs >= {day_obs_min} AND visit1.day_obs <= {day_obs_max}\n", + " -- AND visit1.science_program in ('BLOCK-365', 'BLOCK-407', 'BLOCK-408', 'BLOCK-416', 'BLOCK-417', 'BLOCK-419');\n", + " AND visit1.science_program in ('BLOCK-365', 'BLOCK-407', 'BLOCK-408', 'BLOCK-416', 'BLOCK-417', 'BLOCK-419', 'BLOCK-421', 'BLOCK-T698')\n", + " -- AND visit1.science_program in ('BLOCK-T539', 'BLOCK-T698', 'BLOCK-T641') -- initial alignment, closed loop stability\n", + " AND visit1.physical_filter in ('u_24', 'g_6', 'r_57', 'i_39', 'z_20', 'y_10');\n", + " '''\n", + " \n", + " ccdvisits = client.query(visits_query).to_pandas()\n", + "\n", + " return ccdvisits" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1dd36a19-b1f0-4d2c-81f5-f04faf0c0965", + "metadata": {}, + "outputs": [], + "source": [ + "#day_obs_min = 20251026\n", + "#day_obs_max = 20260106\n", + "zernikes = fetchZernikes(day_obs_min, day_obs_max)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1f400140-a453-48b1-9384-ea5850b390ff", + "metadata": {}, + "outputs": [], + "source": [ + "np.unique(zernikes[\"detector\"])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "75f160c4-0d6c-4950-a458-b5bfc47b664e", + "metadata": {}, + "outputs": [], + "source": [ + "zernike_columns = [f\"z{i}\" for i in range(4, 29)]\n", + "zernikes[\"zernikes\"] = zernikes[zernike_columns].apply(\n", + " lambda row: np.array(row.fillna(0.0).values, dtype=float), axis=1\n", + ")\n", + "zernikes[\"zernikes_fwhm\"] = zernikes[\"zernikes\"].apply(\n", + " convertZernikesToPsfWidth\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "04d4ca0d-d68d-4161-bd49-3ea9053ca485", + "metadata": {}, + "outputs": [], + "source": [ + "zernikes_fwhm = np.vstack(zernikes[\"zernikes_fwhm\"])\n", + "for ii, column in enumerate(zernike_columns):\n", + " zernikes[f\"{column}_fwhm\"] = zernikes_fwhm[:,ii]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "508fee1f-4a09-4794-992c-196e1b835136", + "metadata": {}, + "outputs": [], + "source": [ + "zernikes[\"z4_fwhm\"]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "71f5d9db-483f-434b-b344-dab6a4269aca", + "metadata": {}, + "outputs": [], + "source": [ + "zernikes.columns" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e69e8cf6-7fd9-47fa-9846-f4faf7fcd7ef", + "metadata": {}, + "outputs": [], + "source": [ + "zernikes[['z4', 'detector', 'visit_id']]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "438e3d62-a055-4c19-b668-6d8a731e08dd", + "metadata": {}, + "outputs": [], + "source": [ + "groups = zernikes.groupby('visit_id')\n", + "visits_summary_zernikes = pd.DataFrame({\n", + " 'z4': groups['z4'].mean(),\n", + " 'z4_fwhm': groups['z4_fwhm'].mean(),\n", + "})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "04fa0a07-f11f-4f3a-beda-4c28b76b3c70", + "metadata": {}, + "outputs": [], + "source": [ + "visits_summary_zernikes" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "72987a1d-16a9-49a7-8971-1316104e95d2", + "metadata": {}, + "outputs": [], + "source": [ + "#zernikes[['visit_id', 'z4']][zernikes['detector'] == 195].rename(columns={'z4': 'z4_191'})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7a67cf4c-06aa-43e1-88cc-daa4cf787287", + "metadata": {}, + "outputs": [], + "source": [ + "visits_summary_zernikes = visits_summary_zernikes.merge(zernikes[['visit_id', 'z4']][zernikes['detector'] == 191].rename(columns={'z4': 'z4_191'}), on='visit_id')\n", + "visits_summary_zernikes = visits_summary_zernikes.merge(zernikes[['visit_id', 'z4']][zernikes['detector'] == 195].rename(columns={'z4': 'z4_195'}), on='visit_id')\n", + "visits_summary_zernikes = visits_summary_zernikes.merge(zernikes[['visit_id', 'z4']][zernikes['detector'] == 199].rename(columns={'z4': 'z4_199'}), on='visit_id')\n", + "visits_summary_zernikes = visits_summary_zernikes.merge(zernikes[['visit_id', 'z4']][zernikes['detector'] == 203].rename(columns={'z4': 'z4_203'}), on='visit_id')\n", + "visits_summary_zernikes" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "db3a5e48-7af8-405d-92c9-d927206bfa1d", + "metadata": {}, + "outputs": [], + "source": [ + "\"\"\"\n", + "lists = [np.unique(zernikes['visit_id'][zernikes['detector'] == 191]),\n", + " np.unique(zernikes['visit_id'][zernikes['detector'] == 195]),\n", + " np.unique(zernikes['visit_id'][zernikes['detector'] == 199]),\n", + " np.unique(zernikes['visit_id'][zernikes['detector'] == 203])]\n", + "result = np.array(list(set(lists[0]).intersection(*lists[1:])))\n", + "len(result)\n", + "\"\"\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "32c45459-6c8f-48e1-ab18-639d803a3687", + "metadata": {}, + "outputs": [], + "source": [ + "#visits_summary_zernikes = visits_summary_zernikes[np.isin(visits_summary_zernikes.index, result)]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cd91048b-40c0-449d-a7e1-b2773915fc0d", + "metadata": {}, + "outputs": [], + "source": [ + "\"\"\"\n", + "selection = np.isin(zernikes.index, result)\n", + "visits_summary_zernikes['z4_191'] = zernikes['z4'][zernikes['detector'] == 191].values\n", + "visits_summary_zernikes['z4_195'] = zernikes['z4'][zernikes['detector'] == 195].values\n", + "visits_summary_zernikes['z4_199'] = zernikes['z4'][zernikes['detector'] == 199].values\n", + "visits_summary_zernikes['z4_203'] = zernikes['z4'][zernikes['detector'] == 203].values\n", + "\"\"\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "367e4019-048e-44b2-8ebe-f5593cb4e178", + "metadata": {}, + "outputs": [], + "source": [ + "#visits_summary_zernikes" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4a51d727-0ba7-4392-a612-2b9fe585f410", + "metadata": {}, + "outputs": [], + "source": [ + "#zernikes['z4'][zernikes['detector'] == 191]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "43908b7d-ce54-4951-a67c-69484d6dd48b", + "metadata": {}, + "outputs": [], + "source": [ + "#plt.figure()\n", + "#plt.hist(visits_summary_zernikes['z4_fwhm'], bins=np.linspace(-1., 1., 81))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "38f0a48a-bcb9-486e-836f-f2a78a1cc116", + "metadata": {}, + "outputs": [], + "source": [ + "print(len(visits_summary_zernikes))\n", + "print(len(visits_summary))\n", + "#assert len(visits_summary_zernikes) == len(visits_summary)\n", + "\n", + "visits_summary = pd.merge(visits_summary, visits_summary_zernikes, on=\"visit_id\")\n", + "print(len(visits_summary))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "574d3813-56df-4081-9636-a969a019c29d", + "metadata": {}, + "outputs": [], + "source": [ + "visits_summary.columns" + ] + }, + { + "cell_type": "markdown", + "id": "9c264ca2-3d60-44ee-a834-7167922dc062", + "metadata": { + "jp-MarkdownHeadingCollapsed": true + }, + "source": [ + "# Focus" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1873175a-fac2-4513-87d9-070ea794c567", + "metadata": {}, + "outputs": [], + "source": [ + "bins = np.linspace(-1., 1., 81)\n", + "\n", + "selection = (visits_summary['science_program'] == 'BLOCK-T698')\n", + "\n", + "print(visits_summary['day_obs'][selection].value_counts().sort_index().to_string())\n", + "\n", + "plt.figure()\n", + "#plt.hist(zernikes[\"z4\"], bins=bins, histtype='step')\n", + "#plt.hist(visits_summary[\"z4\"], bins=np.linspace(-1., 1., 41), histtype='step')\n", + "plt.hist(visits_summary[\"z4_191\"][selection], bins=bins, histtype='step', label='191')#, color='tab:blue')\n", + "plt.hist(visits_summary[\"z4_195\"][selection], bins=bins, histtype='step', label='195')\n", + "plt.hist(visits_summary[\"z4_199\"][selection], bins=bins, histtype='step', label='199')\n", + "plt.hist(visits_summary[\"z4_203\"][selection], bins=bins, histtype='step', label='203')\n", + "plt.xlim(-0.5 - 0.25, 0.5 - 0.25)\n", + "plt.legend()\n", + "plt.xlabel('Z4')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "18a0dd2b-0aff-482d-ac00-21d9a9095cd4", + "metadata": {}, + "outputs": [], + "source": [ + "bins = np.linspace(-1., 1., 81)\n", + "selection = (visits_summary['science_program'] == 'BLOCK-T698')\n", + "\n", + "plt.figure()\n", + "for day_obs in np.unique(visits_summary[\"day_obs\"][selection]):\n", + " selection_day_obs = visits_summary[\"day_obs\"] == day_obs\n", + " plt.hist(visits_summary[\"z4_203\"][selection & selection_day_obs], bins=bins, histtype='step', label='203')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0b014fe4-c499-4ccc-9254-d900617768fc", + "metadata": {}, + "outputs": [], + "source": [ + "bins = np.linspace(-1., 1., 81)\n", + "selection = (visits_summary['science_program'] == 'BLOCK-T698')\n", + "\n", + "fig, ax = plt.subplots(4, 1, figsize=(6, 8), sharex=False)\n", + "#fig.subplots_adjust(hspace=None)\n", + "\n", + "for ii, day_obs in enumerate(np.unique(visits_summary[\"day_obs\"][selection])):\n", + " selection_day_obs = visits_summary[\"day_obs\"] == day_obs\n", + " ax[ii].hist(visits_summary[\"z4_191\"][selection & selection_day_obs] - visits_summary[\"z4\"][selection & selection_day_obs], bins=bins, histtype='step', label='191', color='tab:blue')\n", + " ax[ii].hist(visits_summary[\"z4_195\"][selection & selection_day_obs] - visits_summary[\"z4\"][selection & selection_day_obs], bins=bins, histtype='step', label='195', color='tab:orange')\n", + " ax[ii].hist(visits_summary[\"z4_199\"][selection & selection_day_obs] - visits_summary[\"z4\"][selection & selection_day_obs], bins=bins, histtype='step', label='199', color='tab:green')\n", + " ax[ii].hist(visits_summary[\"z4_203\"][selection & selection_day_obs] - visits_summary[\"z4\"][selection & selection_day_obs], bins=bins, histtype='step', label='203', color='tab:red')\n", + " ax[ii].set_xlim(-0.5, 0.5)\n", + " ax[ii].set_xlabel(r'Z4_sensor - Z4_average')\n", + " ax[ii].set_ylabel('Visit Counts')\n", + " #ax[ii].set_title(str(day_obs))\n", + " ax[ii].legend(title=str(day_obs))\n", + "\n", + "plt.tight_layout()\n", + "plt.show()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cf0a8914-ccac-4903-8c30-0cb271585af0", + "metadata": {}, + "outputs": [], + "source": [ + "print(np.median(visits_summary[\"z4_191\"]))\n", + "print(np.median(visits_summary[\"z4_195\"]))\n", + "print(np.median(visits_summary[\"z4_199\"]))\n", + "print(np.median(visits_summary[\"z4_203\"]))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9d577b30-4874-4a28-9a2f-db8125b34ffb", + "metadata": {}, + "outputs": [], + "source": [ + "#x = np.std(np.vstack([\n", + "# visits_summary[\"z4_191\"] - np.median(visits_summary[\"z4_191\"]),\n", + "# visits_summary[\"z4_195\"] - np.median(visits_summary[\"z4_195\"]),\n", + "# visits_summary[\"z4_199\"] - np.median(visits_summary[\"z4_199\"]),\n", + "# visits_summary[\"z4_203\"] - np.median(visits_summary[\"z4_203\"]),\n", + "#]), axis=0)\n", + "x = np.std(np.vstack([\n", + " visits_summary[\"z4_191\"],\n", + " visits_summary[\"z4_195\"],\n", + " visits_summary[\"z4_199\"],\n", + " visits_summary[\"z4_203\"],\n", + "]), axis=0)\n", + "plt.figure()\n", + "#plt.hist(x, bins=np.linspace(0., 0.4, 41))\n", + "plt.scatter(x, visits_summary[\"z4\"], s=1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e0440371-0e88-4a76-b785-f58dc43e0dd4", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "plt.scatter(x, visits_summary[\"psf_fwhm_95_05\"], s=1, c= visits_summary['psf_fwhm_50'], vmin=0.6, vmax=1.2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0e0b6a01-95eb-43d0-a7dc-c08dc960d0a6", + "metadata": {}, + "outputs": [], + "source": [ + "bins = np.linspace(-1., 1., 81)\n", + "\n", + "selection_prelsst = (visits_summary['science_program'] == 'BLOCK-407')\n", + "selection_stability = (visits_summary['science_program'] == 'BLOCK-T698')\n", + "\n", + "plt.figure()\n", + "#plt.hist(zernikes[\"z4\"], bins=bins, histtype='step')\n", + "#plt.hist(visits_summary[\"z4\"], bins=np.linspace(-1., 1., 41), histtype='step')\n", + "plt.hist(visits_summary[\"z4\"][selection_prelsst], bins=bins, histtype='step', density=True, label='BLOCK-407')\n", + "plt.hist(visits_summary[\"z4\"][selection_stability], bins=bins, histtype='step', density=True, label='BLOCK-T698')\n", + "plt.legend()\n", + "plt.xlim(-1., 1.)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c5f20dae-9168-4c90-a9ca-23e4ce82f246", + "metadata": {}, + "outputs": [], + "source": [ + "day_obs = 20260408\n", + "selection = (visits_summary['day_obs'] == day_obs)\n", + "\n", + "plt.figure(figsize=(16, 4))\n", + "plt.axhline(-0.3, ls='--', c='0.8')\n", + "plt.axhline(+0.3, ls='--', c='0.8')\n", + "#plt.scatter(visits_summary[\"seq_num\"], visits_summary[\"z4_191\"])\n", + "#plt.scatter(visits_summary[\"seq_num\"], visits_summary[\"z4_195\"])\n", + "#plt.scatter(visits_summary[\"seq_num\"], visits_summary[\"z4_199\"])\n", + "#plt.scatter(visits_summary[\"seq_num\"], visits_summary[\"z4_203\"])\n", + "#plt.plot(visits_summary[\"seq_num\"][selection], visits_summary[\"z4_191\"][selection])\n", + "#plt.plot(visits_summary[\"seq_num\"][selection], visits_summary[\"z4_195\"][selection])\n", + "#plt.plot(visits_summary[\"seq_num\"][selection], visits_summary[\"z4_199\"][selection])\n", + "#plt.plot(visits_summary[\"seq_num\"][selection], visits_summary[\"z4_203\"][selection])\n", + "#plt.plot(visits_summary[\"seq_num\"][selection], -0.5 + visits_summary[\"z4\"][selection], label='Z4_mean', c='black')\n", + "plt.plot(visits_summary[\"seq_num\"][selection], +0.3 + (visits_summary[\"z4_191\"][selection] - visits_summary[\"z4\"][selection]), label='191 (offset by +0.3 for visualization)')\n", + "plt.plot(visits_summary[\"seq_num\"][selection], -0.3 + (visits_summary[\"z4_195\"][selection] - visits_summary[\"z4\"][selection]), label='195 (offset by -0.3 for visualization)')\n", + "plt.plot(visits_summary[\"seq_num\"][selection], -0.3 + (visits_summary[\"z4_199\"][selection] - visits_summary[\"z4\"][selection]), label='199 (offset by -0.3 for visualization)')\n", + "plt.plot(visits_summary[\"seq_num\"][selection], +0.3 + (visits_summary[\"z4_203\"][selection] - visits_summary[\"z4\"][selection]), label='203 (offset by +0.3 for visualization)')\n", + "plt.legend()\n", + "plt.xlabel('seq_num')\n", + "plt.ylabel('Z4_sensor - Z4_mean')\n", + "plt.title(day_obs)\n", + "\n", + "#ax2 = plt.gca().twinx()\n", + "#ax2.plot(visits_summary[\"seq_num\"][selection], -0.5 + visits_summary[\"z4\"][selection], label='Z4_mean', c='black')\n", + "#ax2.set_label('Z4_mean')\n", + "#plt.ylim(-0.75, 0.25)\n", + "\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7061cd40-0d2e-4a57-9ccb-3310e643bb7e", + "metadata": {}, + "outputs": [], + "source": [ + "selection_prelsst = (visits_summary['science_program'] == 'BLOCK-407')\n", + "selection_stability = (visits_summary['science_program'] == 'BLOCK-T698')\n", + "\n", + "bins = np.linspace(-0.4, 0.4, 21)\n", + "\n", + "values_prelsst = np.concatenate([\n", + " visits_summary[\"z4_191\"][selection & selection_prelsst] - visits_summary[\"z4\"][selection & selection_prelsst],\n", + " visits_summary[\"z4_195\"][selection & selection_prelsst] - visits_summary[\"z4\"][selection & selection_prelsst],\n", + " visits_summary[\"z4_199\"][selection & selection_prelsst] - visits_summary[\"z4\"][selection & selection_prelsst],\n", + " visits_summary[\"z4_203\"][selection & selection_prelsst] - visits_summary[\"z4\"][selection & selection_prelsst],\n", + "])\n", + "values_stability = np.concatenate([\n", + " visits_summary[\"z4_191\"][selection & selection_stability] - visits_summary[\"z4\"][selection & selection_stability],\n", + " visits_summary[\"z4_195\"][selection & selection_stability] - visits_summary[\"z4\"][selection & selection_stability],\n", + " visits_summary[\"z4_199\"][selection & selection_stability] - visits_summary[\"z4\"][selection & selection_stability],\n", + " visits_summary[\"z4_203\"][selection & selection_stability] - visits_summary[\"z4\"][selection & selection_stability],\n", + "])\n", + "\n", + "plt.figure()\n", + "plt.hist(visits_summary[\"z4_199\"][selection & selection_prelsst] - visits_summary[\"z4\"][selection & selection_prelsst], bins=bins, histtype='step')\n", + "plt.hist(visits_summary[\"z4_199\"][selection & selection_stability] - visits_summary[\"z4\"][selection & selection_stability], bins=bins, histtype='step')\n", + "#plt.hist(values_prelsst, bins=bins, histtype='step')\n", + "#plt.hist(values_stability, bins=bins, histtype='step')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0f98db35-6b64-478a-a0c9-bdf12c35c5e9", + "metadata": {}, + "outputs": [], + "source": [ + "np.concatenate([\n", + " [0, 1, 2],\n", + " [3, 4, 5],\n", + "])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9c5ade07-5066-426c-ab04-b48fb7bcc6aa", + "metadata": {}, + "outputs": [], + "source": [ + "selection = (visits_summary['donut_blur_fwhm'] < 0.9)\n", + "\n", + "plt.figure()\n", + "plt.scatter(#np.fabs(visits_summary[\"z4_fwhm\"] + 0.25), \n", + " #visits_summary[\"z4_fwhm\"],\n", + " #visits_summary[\"aos_fwhm\"],\n", + " visits_summary[\"z4\"][selection],\n", + " visits_summary[\"psf_fwhm_50\"][selection], \n", + " #visits_summary[\"psf_fwhm_95_05\"],\n", + " #visits_summary[\"psf_e_50\"], \n", + " s=4, \n", + " c=visits_summary['donut_blur_fwhm'][selection], vmin=0.7, vmax=1.2, cmap='viridis')\n", + "plt.colorbar(label='Donut Blur')\n", + "plt.xlabel('Mean z4')\n", + "plt.ylabel('Median PSF FWHM (arcsec)')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "75b66450-0dd5-4375-939e-7a45a70f4b05", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "plt.scatter(visits_summary[\"z4_fwhm\"],\n", + " visits_summary[\"aos_fwhm\"],\n", + " s=4) \n", + " #c=visits_summary['psf_e_50'], vmin=0.05, vmax=0.2)\n", + "#plt.colorbar(label='Donut Blur')\n", + "plt.xlabel('Mean z4')\n", + "plt.ylabel('aos_fwhm')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cd5b1666-2bf2-448d-b25d-03e5f384754d", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "plt.hist(visits_summary[\"z4\"], bins=)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "eace8b7d-2a0a-4a8b-b801-bc8cb9d0c130", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "c04288cf-f976-488c-b790-0c8e960fd782", + "metadata": { + "jp-MarkdownHeadingCollapsed": true + }, + "source": [ + "## Correlation Analysis" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "92dd65c1-23f0-4336-8f18-ec63273e1b57", + "metadata": {}, + "outputs": [], + "source": [ + "#plt.figure()\n", + "#plt.scatter(visits_summary[\"physical_rotator_angle\"], visits_summary[\"psf_fwhm_95_05\"], s=2, edgecolor='none')\n", + "#plt.scatter(visits_summary[\"physical_rotator_angle\"], visits_summary[\"aos_fwhm\"], s=2, edgecolor='none')\n", + "#plt.scatter(visits_summary[\"physical_rotator_angle\"], visits_summary[\"sys_50\"], s=2, edgecolor='none')\n", + "#plt.scatter(visits_summary[\"physical_rotator_angle\"], visits_summary[\"psf_fwhm_50\"], s=2, edgecolor='none')\n", + "#plt.hexbin(visits_summary[\"physical_rotator_angle\"], visits_summary[\"psf_fwhm_95_05\"], gridsize=20, mincnt=1,)\n", + "#plt.hexbin(visits_summary[\"physical_rotator_angle\"], visits_summary[\"aos_fwhm\"], gridsize=20, mincnt=1,)\n", + "#plt.hexbin(visits_summary[\"physical_rotator_angle\"], visits_summary[\"sys_50\"], gridsize=20, mincnt=1,)\n", + "\n", + "metrics = ('aos_cam_fwhm', 'sys_50', 'psf_fwhm_95_05')\n", + "fig, ax = plt.subplots(len(metrics), sharex=True, figsize=(6, 8))\n", + "\n", + "for ii, metric in enumerate(metrics):\n", + " print(ii, metric)\n", + " hb = ax[ii].hexbin(visits_summary[\"physical_rotator_angle\"], visits_summary[metric], gridsize=20, mincnt=1,)\n", + " ax[ii].set_xlim(-90., 90.)\n", + " ax[ii].set_ylim(0., 2.)\n", + " #ax[ii].set_title(metric)\n", + " ax[ii].set_xlabel('Phys Rot Angle (deg)')\n", + " ax[ii].set_ylabel('%s (arcsec)'%(metric))\n", + " fig.colorbar(hb, ax=ax[ii], label='Counts')\n", + " \n", + "#ax[0].legend()\n", + "#ax[-1].set_xlabel('%s (arsec)'%(metric))\n", + "#plt.suptitle('Physical Rotator Angle Changes')\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f1eb047e-e303-4d37-bb52-1fc00d9596b7", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "plt.scatter(visits_summary[\"physical_rotator_angle\"], visits_summary[\"altitude\"], c=visits_summary[\"psf_fwhm_95_05\"], vmin=0.3, vmax=0.6, s=5)\n", + "plt.colorbar()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "769fd484-bfa3-4eef-9eae-8db836271cf9", + "metadata": {}, + "outputs": [], + "source": [ + "#plt.figure()\n", + "#plt.scatter(visits_summary[\"altitude\"], visits_summary[\"psf_fwhm_95_05\"], s=2, edgecolor='none')\n", + "#plt.scatter(visits_summary[\"altitude\"], visits_summary[\"aos_fwhm\"], s=2, edgecolor='none')\n", + "#plt.scatter(visits_summary[\"altitude\"], visits_summary[\"sys_50\"], s=2, edgecolor='none')\n", + "#plt.scatter(visits_summary[\"altitude\"], visits_summary[\"psf_fwhm_50\"], s=2, edgecolor='none')\n", + "#plt.hexbin(visits_summary[\"altitude\"], visits_summary[\"psf_fwhm_95_05\"], gridsize=20, mincnt=1,)\n", + "#plt.hexbin(visits_summary[\"altitude\"], visits_summary[\"aos_fwhm\"], gridsize=20, mincnt=1,)\n", + "#plt.hexbin(visits_summary[\"altitude\"], visits_summary[\"sys_50\"], gridsize=20, mincnt=1,)\n", + "\n", + "metrics = ('aos_cam_fwhm', 'sys_50', 'psf_fwhm_95_05')\n", + "fig, ax = plt.subplots(len(metrics), sharex=True, figsize=(6, 8))\n", + "\n", + "for ii, metric in enumerate(metrics):\n", + " print(ii, metric)\n", + " hb = ax[ii].hexbin(visits_summary[\"altitude\"], visits_summary[metric], gridsize=20, mincnt=1,)\n", + " ax[ii].set_xlim(30., 90.)\n", + " ax[ii].set_ylim(0., 2.)\n", + " #ax[ii].set_title(metric)\n", + " ax[ii].set_xlabel('Altitude (deg)')\n", + " ax[ii].set_ylabel('%s (arcsec)'%(metric))\n", + " fig.colorbar(hb, ax=ax[ii], label='Counts')\n", + " \n", + "#ax[0].legend()\n", + "#ax[-1].set_xlabel('%s (arsec)'%(metric))\n", + "#plt.suptitle('Physical Rotator Angle Changes')\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3011ee4a-1789-41ca-a8a5-8cc5f7ddc56b", + "metadata": {}, + "outputs": [], + "source": [ + "\"\"\"\n", + "selection = (visits_summary[\"altitude\"] > 40.)\n", + "\n", + "bins = np.linspace(0.2, 2, 21)\n", + "\n", + "#metric = 'aos_fwhm'\n", + "metric = 'sys_50'\n", + "#metric = 'psf_fwhm_95_05'\n", + "#metric = 'psf_fwhm_50'\n", + "\n", + "plt.figure()\n", + "plt.hist(visits_summary[metric][selection], bins=bins, density=True, histtype='step')\n", + "plt.hist(visits_summary[metric][~selection], bins=bins, density=True, histtype='step')\n", + "\"\"\"\n", + "\n", + "\n", + "threshold = 40\n", + "selection = (visits_summary[\"altitude\"] < threshold)\n", + "bins = np.linspace(0.2, 2., 21)\n", + "\n", + "#metric = 'aos_fwhm'\n", + "#metric = 'sys_50'\n", + "#metric = 'psf_fwhm_95_05'\n", + "#metric = 'psf_fwhm_50'\n", + "\n", + "#plt.figure()\n", + "#plt.hist(visits_summary[metric][1:][selection], bins=bins, histtype='step', density=True, label='Filter Change (n = %i)'%(np.sum(selection)))\n", + "#plt.hist(visits_summary[metric][1:][~selection], bins=bins, histtype='step', density=True, label='No Filter Change (n = %i)'%(np.sum(~selection)))\n", + "#plt.legend()\n", + "\n", + "metrics = ('aos_cam_fwhm', 'sys_50', 'psf_fwhm_95_05')\n", + "fig, ax = plt.subplots(len(metrics), sharex=False, figsize=(6, 8))\n", + "\n", + "for ii, metric in enumerate(metrics):\n", + " print(ii, metric)\n", + " ax[ii].hist(visits_summary[metric][1:][selection], bins=bins, histtype='step', density=True, label='Altitude < %.f deg (n = %i)'%(threshold, np.sum(selection)))\n", + " ax[ii].hist(visits_summary[metric][1:][~selection], bins=bins, histtype='step', density=True, label='Altitude > %.f deg (n = %i)'%(threshold, np.sum(~selection)))\n", + " ax[ii].set_xlim(0.2, 2.)\n", + " #ax[ii].set_title(metric)\n", + " ax[ii].set_xlabel('%s (arsec)'%(metric))\n", + " ax[ii].set_ylabel('PDF')\n", + " \n", + "ax[0].legend()\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c858006e-d3f4-4282-b793-71c16ce5b1f8", + "metadata": {}, + "outputs": [], + "source": [ + "threshold = 0.\n", + "selection = (visits_summary[\"physical_rotator_angle\"] < threshold)\n", + "bins = np.linspace(0.2, 2., 21)\n", + "\n", + "metrics = ('aos_cam_fwhm', 'sys_50', 'psf_fwhm_95_05')\n", + "fig, ax = plt.subplots(len(metrics), sharex=False, figsize=(6, 8))\n", + "\n", + "for ii, metric in enumerate(metrics):\n", + " print(ii, metric)\n", + " ax[ii].hist(visits_summary[metric][1:][selection], bins=bins, histtype='step', density=True, label='Phys Rot Angle < %.f deg (n = %i)'%(threshold, np.sum(selection)))\n", + " ax[ii].hist(visits_summary[metric][1:][~selection], bins=bins, histtype='step', density=True, label='Phys Rot Angle > %.f deg (n = %i)'%(threshold, np.sum(~selection)))\n", + " ax[ii].set_xlim(0.2, 2.)\n", + " #ax[ii].set_title(metric)\n", + " ax[ii].set_xlabel('%s (arsec)'%(metric))\n", + " ax[ii].set_ylabel('PDF')\n", + " \n", + "ax[0].legend()\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "71862e1e-02d6-49d9-8704-973525f11f96", + "metadata": {}, + "outputs": [], + "source": [ + "selection = (visits_summary['day_obs'] < 20251126)\n", + "\n", + "bins = np.linspace(20, 90, 21)\n", + "\n", + "plt.figure()\n", + "plt.hist(visits_summary['altitude'][selection], bins=bins, density=True, histtype='step')\n", + "plt.hist(visits_summary['altitude'][~selection], bins=bins, density=True, histtype='step')" + ] + }, + { + "cell_type": "markdown", + "id": "c2e12ab1-b90f-4228-9f22-d7bb9d5f233f", + "metadata": { + "jp-MarkdownHeadingCollapsed": true + }, + "source": [ + "## DDF" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "183ace4e-e585-4ddb-9929-3946deeb3455", + "metadata": {}, + "outputs": [], + "source": [ + "selection = visits_summary['target_name'].str.contains('ddf')\n", + "\n", + "plt.figure()\n", + "for band in bands:\n", + " band_selection = selection & (visits_summary[\"band\"] == band)\n", + " plt.scatter(\n", + " visits_summary[\"exp_midpt_mjd\"][selection & band_selection], \n", + " visits_summary[\"aos_fwhm\"][selection & band_selection],\n", + " c=colors[band], label=band,\n", + " edgecolor='none', s=10\n", + " )\n", + "plt.legend()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3974b311-64c1-4b59-8b00-78c9a1483ff9", + "metadata": {}, + "outputs": [], + "source": [ + "selection = visits_summary['target_name'].str.contains('ddf')\n", + "\n", + "plt.figure()\n", + "for band in bands:\n", + " band_selection = selection & (visits_summary[\"band\"] == band)\n", + " plt.scatter(visits_summary['altitude'][selection & band_selection], visits_summary[''][selection & band_selection], c=colors[band], label=band,)\n", + "plt.legend()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bf51e04b-2f4c-47c4-a3e2-0bf95d136c71", + "metadata": {}, + "outputs": [], + "source": [ + "visits_summary[selection]['day_obs'].value_counts().sort_index()" + ] + }, + { + "cell_type": "markdown", + "id": "8f42d4a0-6399-4294-a853-a8cf22d6cf26", + "metadata": { + "jp-MarkdownHeadingCollapsed": true + }, + "source": [ + "## Differences" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2264af1f-2bc8-4eeb-93b2-55fd7c649fa8", + "metadata": {}, + "outputs": [], + "source": [ + "a = [1, 10, 100, 1000]\n", + "np.diff(a)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b859c727-9d2e-4a15-9f84-3d8c58ed5a07", + "metadata": {}, + "outputs": [], + "source": [ + "selection_filter_change = (visits_summary[\"band\"][1:].values != visits_summary[\"band\"][0:-1].values)\n", + "selection_consecutive = np.diff(visits_summary['seq_num']) == 1\n", + "selection = (selection_consecutive & selection_filter_change)\n", + "\n", + "bins = np.linspace(0.2, 2.0, 41)\n", + "\n", + "#metric = 'aos_fwhm'\n", + "#metric = 'sys_50'\n", + "#metric = 'psf_fwhm_95_05'\n", + "#metric = 'psf_fwhm_50'\n", + "\n", + "#plt.figure()\n", + "#plt.hist(visits_summary[metric][1:][selection], bins=bins, histtype='step', density=True, label='Filter Change (n = %i)'%(np.sum(selection)))\n", + "#plt.hist(visits_summary[metric][1:][~selection], bins=bins, histtype='step', density=True, label='No Filter Change (n = %i)'%(np.sum(~selection)))\n", + "#plt.legend()\n", + "\n", + "metrics = ('aos_cam_fwhm', 'sys_50', 'psf_fwhm_95_05')\n", + "fig, ax = plt.subplots(len(metrics), sharex=False, figsize=(6, 8))\n", + "\n", + "for ii, metric in enumerate(metrics):\n", + " print(ii, metric)\n", + " ax[ii].hist(visits_summary[metric][1:][selection], bins=bins, histtype='step', density=True, label='Filter Change (n = %i)'%(np.sum(selection)))\n", + " ax[ii].hist(visits_summary[metric][1:][~selection], bins=bins, histtype='step', density=True, label='No Filter Change (n = %i)'%(np.sum(~selection)))\n", + " ax[ii].set_xlim(0.2, 2.)\n", + " #ax[ii].set_title(metric)\n", + " ax[ii].set_xlabel('%s (arsec)'%(metric))\n", + " ax[ii].set_ylabel('PDF')\n", + " \n", + "ax[0].legend()\n", + "#ax[-1].set_xlabel('%s (arsec)'%(metric))\n", + "#plt.suptitle('Physical Rotator Angle Changes')\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c9aeab7c-6ffb-4aa0-95e3-227647073309", + "metadata": {}, + "outputs": [], + "source": [ + "\"\"\"\n", + "selection_consecutive = np.diff(visits_summary['seq_num']) == 1\n", + "\n", + "#metric = 'aos_fwhm'\n", + "#metric = 'sys_50'\n", + "metric = 'psf_fwhm_95_05'\n", + "#metric = 'psf_fwhm_50'\n", + "\n", + "plt.figure()\n", + "plt.scatter(np.diff(visits_summary[\"altitude\"])[selection_consecutive], visits_summary[metric][1:][selection_consecutive], s=1)\n", + "#plt.scatter(np.diff(visits_summary[\"physical_rotator_angle\"])[selection_consecutive], visits_summary[metric][1:][selection_consecutive], s=1)\n", + "\"\"\"\n", + "\n", + "threshold = 10.\n", + "\n", + "metrics = ('aos_cam_fwhm', 'sys_50', 'psf_fwhm_95_05')\n", + "fig, ax = plt.subplots(len(metrics), sharex=False, figsize=(6, 8))\n", + "\n", + "for ii, metric in enumerate(metrics):\n", + " print(ii, metric)\n", + " ax[ii].axvspan(-1. * threshold, threshold, color='0.9')\n", + " ax[ii].scatter(np.diff(visits_summary[\"altitude\"])[selection_consecutive], visits_summary[metric][1:][selection_consecutive], s=2, edgecolor='none')\n", + " ax[ii].set_xlim(-45, 45.)\n", + " ax[ii].set_ylim(0.2, 2.)\n", + " #ax[ii].set_title(metric)\n", + " ax[ii].set_xlabel(r'$\\Delta$ altitude (deg)')\n", + " ax[ii].set_ylabel('%s (arsec)'%(metric))\n", + "\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "079a5ccc-9ada-4af2-949a-85c7ed552e36", + "metadata": {}, + "outputs": [], + "source": [ + "threshold = 10. \n", + "\n", + "metrics = ('aos_cam_fwhm', 'sys_50', 'psf_fwhm_95_05')\n", + "fig, ax = plt.subplots(len(metrics), sharex=False, figsize=(6, 8))\n", + "\n", + "for ii, metric in enumerate(metrics):\n", + " print(ii, metric)\n", + " ax[ii].axvspan(-1. * threshold, threshold, color='0.9')\n", + " ax[ii].scatter(np.diff(visits_summary[\"physical_rotator_angle\"])[selection_consecutive], visits_summary[metric][1:][selection_consecutive], s=2, edgecolor='none')\n", + " ax[ii].set_xlim(-100, 100.)\n", + " ax[ii].set_ylim(0.2, 2.)\n", + " #ax[ii].set_title(metric)\n", + " ax[ii].set_xlabel(r'$\\Delta$ phys rot angle (deg)')\n", + " ax[ii].set_ylabel('%s (arsec)'%(metric))\n", + "\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7e9c17ee-4988-4db7-b4a9-3c600c00e827", + "metadata": {}, + "outputs": [], + "source": [ + "selection_consecutive = np.diff(visits_summary['seq_num']) == 1\n", + "#selection_slew = (np.fabs(np.diff(visits_summary[\"altitude\"])) > 10)\n", + "selection_slew = (np.fabs(np.diff(visits_summary[\"physical_rotator_angle\"])) > 20.)\n", + "\n", + "selection = (selection_slew & selection_consecutive)\n", + "\n", + "bins = np.linspace(0.2, 2.0, 41)\n", + "\n", + "metric = 'aos_fwhm'\n", + "#metric = 'sys_50'\n", + "#metric = 'psf_fwhm_95_05'\n", + "#metric = 'psf_fwhm_50'\n", + "\n", + "plt.figure()\n", + "plt.hist(visits_summary[metric][1:][selection], bins=bins, histtype='step', density=True, label='Large Slew (n = %i)'%(np.sum(selection)))\n", + "plt.hist(visits_summary[metric][1:][~selection], bins=bins, histtype='step', density=True, label='No large Slew (n = %i)'%(np.sum(~selection)))\n", + "plt.legend()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f4382141-e4b5-48d7-b6f6-bf4cbe435d54", + "metadata": {}, + "outputs": [], + "source": [ + "selection_consecutive = np.diff(visits_summary['seq_num']) == 1\n", + "#selection_slew = (np.fabs(np.diff(visits_summary[\"altitude\"])) > 10)\n", + "threshold = 10.\n", + "selection_slew = (np.fabs(np.diff(visits_summary[\"physical_rotator_angle\"])) > threshold)\n", + "\n", + "selection = (selection_slew & selection_consecutive)\n", + "\n", + "bins = np.linspace(0.2, 2.0, 41)\n", + "\n", + "metrics = ('aos_cam_fwhm', 'sys_50', 'psf_fwhm_95_05')\n", + "fig, ax = plt.subplots(len(metrics), sharex=False, figsize=(6, 8))\n", + "\n", + "for ii, metric in enumerate(metrics):\n", + " print(ii, metric)\n", + " ax[ii].hist(visits_summary[metric][1:][selection], bins=bins, histtype='step', density=True, label=r'$|\\Delta~\\mathrm{phys~rot~angle}| > %.f$ deg (n = %i)'%(threshold , np.sum(selection)))\n", + " ax[ii].hist(visits_summary[metric][1:][~selection], bins=bins, histtype='step', density=True, label=r'$|\\Delta~\\mathrm{phys~rot~angle}| < %.f$ deg (n = %i)'%(threshold, np.sum(~selection)))\n", + " ax[ii].set_xlim(0.2, 2.)\n", + " #ax[ii].set_title(metric)\n", + " ax[ii].set_xlabel('%s (arsec)'%(metric))\n", + " ax[ii].set_ylabel('PDF')\n", + " \n", + "ax[0].legend()\n", + "#ax[-1].set_xlabel('%s (arsec)'%(metric))\n", + "#plt.suptitle('Physical Rotator Angle Changes')\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "27c2b116-a19a-4829-847a-a2fdd8b7f1f0", + "metadata": {}, + "outputs": [], + "source": [ + "selection_consecutive = np.diff(visits_summary['seq_num']) == 1\n", + "threshold = 10\n", + "selection_slew = (np.fabs(np.diff(visits_summary[\"altitude\"])) > threshold)\n", + "#threshold = 20.\n", + "#selection_slew = (np.fabs(np.diff(visits_summary[\"physical_rotator_angle\"])) > threshold)\n", + "\n", + "selection = (selection_slew & selection_consecutive)\n", + "\n", + "bins = np.linspace(0.2, 2.0, 41)\n", + "\n", + "metrics = ('aos_cam_fwhm', 'sys_50', 'psf_fwhm_95_05')\n", + "fig, ax = plt.subplots(len(metrics), sharex=False, figsize=(6, 8))\n", + "\n", + "for ii, metric in enumerate(metrics):\n", + " print(ii, metric)\n", + " ax[ii].hist(visits_summary[metric][1:][selection], bins=bins, histtype='step', density=True, label=r'$|\\Delta~\\mathrm{altitude}| > %.f$ deg (n = %i)'%(threshold , np.sum(selection)))\n", + " ax[ii].hist(visits_summary[metric][1:][~selection], bins=bins, histtype='step', density=True, label=r'$|\\Delta~\\mathrm{altitude}| < %.f$ deg (n = %i)'%(threshold, np.sum(~selection)))\n", + " ax[ii].set_xlim(0.2, 2.)\n", + " #ax[ii].set_title(metric)\n", + " ax[ii].set_xlabel('%s (arsec)'%(metric))\n", + "\n", + "ax[0].legend()\n", + "#ax[-1].set_xlabel('%s (arsec)'%(metric))\n", + "#plt.suptitle('Physical Rotator Angle Changes')\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "59ddfc4e-55fa-4120-83a4-4c48a78a1207", + "metadata": { + "jp-MarkdownHeadingCollapsed": true + }, + "source": [ + "## Inserted Visits" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "98c4ac0f-2e02-4624-b4c5-35c3a070893a", + "metadata": {}, + "outputs": [], + "source": [ + "np.sum(visits_summary['observation_reason'].str.contains('filter_change_close_loop'))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bfa6cac8-829f-4cfc-8bb1-af30a7e14e5b", + "metadata": {}, + "outputs": [], + "source": [ + "import scipy" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "386a7f71-398a-49d4-8299-dc3e3bb339f1", + "metadata": {}, + "outputs": [], + "source": [ + "# Identify the first visit of filter_change_close_loop sequence\n", + "condition = visits_summary['observation_reason'].str.contains('filter_change_close_loop') #& visits_summary['science_program'] == 'BLOCK-419'\n", + "\n", + "label, num_features = scipy.ndimage.label(condition)\n", + "index_filter_change_close_loop = np.array([label.tolist().index(_) for _ in range(1, num_features)])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3f1868f0-f471-452f-aada-c5e1e4db26e3", + "metadata": {}, + "outputs": [], + "source": [ + "index_filter_change_close_loop" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a092cd6f-7aa0-4105-91b0-ece4b0d124a2", + "metadata": {}, + "outputs": [], + "source": [ + "visits_summary['observation_reason'].iloc[index_filter_change_close_loop]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b66cf7c1-0698-494e-b2a6-4d560806e006", + "metadata": {}, + "outputs": [], + "source": [ + "visits_summary['science_program'].iloc[index_filter_change_close_loop]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1b368795-17c3-4a5e-acce-e3613166913d", + "metadata": {}, + "outputs": [], + "source": [ + "visits_summary['band'].iloc[index_filter_change_close_loop]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "92010639-e775-4b18-9daf-75b1ef9d3b8c", + "metadata": {}, + "outputs": [], + "source": [ + "visits_summary['band'].iloc[index_filter_change_close_loop - 1]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "eb131386-b284-462b-ba41-6c7e4ee1bd58", + "metadata": {}, + "outputs": [], + "source": [ + "print(len(index_filter_change_close_loop))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b9f5e231-dc5b-4a30-b7aa-1ae45688838b", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "\n", + "nvisits = 16\n", + "\n", + "metric = 'aos_fwhm'\n", + "#metric = 'z4_fwhm'\n", + "#metric = 'psf_fwhm_50'\n", + "#metric = 'psf_fwhm_zenith_500nm_50'\n", + "\n", + "markers = ['o', 'v', '^', 's', 'p', 'h', 'H', '+', 'x', 'd']\n", + "\n", + "x_array = []\n", + "y_array = []\n", + "\n", + "exclude = 'y'\n", + "include = None\n", + "assert (exclude is None or include is None)\n", + "\n", + "for ii in range(0, len(index_filter_change_close_loop)):\n", + " #print(ii)\n", + " #print(visits_summary['target_name'].iloc[index_filter_change_close_loop[ii]])\n", + " #color = 'red' if 'ddf' in visits_summary['target_name'].iloc[index_filter_change_close_loop[ii]] else 'black'\n", + "\n", + " before_filter = visits_summary['band'].iloc[index_filter_change_close_loop[ii] - 1]\n", + " after_filter = visits_summary['band'].iloc[index_filter_change_close_loop[ii]]\n", + " #print(before_filter, after_filter)\n", + "\n", + " \"\"\"\n", + " if before_filter == \"y\" or after_filter == \"y\":\n", + " color = 'red'\n", + " #continue\n", + " else:\n", + " color = 'black'\n", + " continue\n", + " \"\"\"\n", + "\n", + " if exclude is not None:\n", + " if before_filter == exclude or after_filter == exclude:\n", + " continue\n", + " if include is not None:\n", + " if not (before_filter == include or after_filter == include):\n", + " continue\n", + "\n", + " print(before_filter, after_filter)\n", + " \n", + " #print(visits_summary['target_name'].iloc[index_filter_change_close_loop[ii] + np.arange(nvisits)] == visits_summary['target_name'].iloc[index_filter_change_close_loop[ii]])\n", + " delta_altitude = np.fabs(\n", + " visits_summary['altitude'].iloc[index_filter_change_close_loop[ii] + np.arange(nvisits)] - visits_summary['altitude'].iloc[index_filter_change_close_loop[ii]]\n", + " )\n", + " delta_physical_rotator_angle = np.fabs(\n", + " visits_summary['physical_rotator_angle'].iloc[index_filter_change_close_loop[ii] + np.arange(nvisits)] - visits_summary['physical_rotator_angle'].iloc[index_filter_change_close_loop[ii]]\n", + " )\n", + " same_filter = (visits_summary['band'].iloc[index_filter_change_close_loop[ii] + np.arange(nvisits)] == visits_summary['band'].iloc[index_filter_change_close_loop[ii]])\n", + " selection = (delta_altitude < 10.) & (delta_physical_rotator_angle < 10.) & (same_filter)\n", + "\n", + " #print(same_filter[selection])\n", + " #print(delta_altitude[selection])\n", + " #print(delta_physical_rotator_angle[selection])\n", + " #print(visits_summary['band'].iloc[index_filter_change_close_loop[ii] + np.arange(nvisits)])\n", + " #print(same_filter)\n", + " #print(delta_altitude < 10.)\n", + " #print(delta_physical_rotator_angle < 10.)\n", + " \n", + " x = np.arange(nvisits) + np.random.uniform(-0.1, 0.1)\n", + " y = visits_summary[metric].iloc[index_filter_change_close_loop[ii] + np.arange(nvisits)]\n", + " #y = visits_summary[metric].iloc[index_filter_change_close_loop[ii] + np.arange(nvisits)] - visits_summary[metric].iloc[index_filter_change_close_loop[ii] - 1]\n", + " #plt.scatter(x[selection], y[selection], marker=markers[ii % len(markers)])\n", + " plt.plot(x[selection], y[selection], lw=1)\n", + " #plt.plot(x[selection], y[selection], lw=1, c=color)\n", + "\n", + "\n", + " x_array.append(x)\n", + " y_array.append(y.values)\n", + "\n", + "\n", + "x_median = []\n", + "y_median = []\n", + "for ii in range(0, nvisits):\n", + " y = []\n", + " for jj in range(0, len(y_array)):\n", + " if len(y_array[jj]) >= ii:\n", + " y.append(y_array[jj][ii])\n", + "\n", + " #print(ii)\n", + " #print(y)\n", + " \n", + " if len(x) >= 3:\n", + " x_median.append(ii)\n", + " y_median.append(np.nanmedian(y))\n", + " \n", + "\n", + "#print(x_median)\n", + "#print(y_median)\n", + "\n", + "plt.plot(x_median, y_median, c='black', lw=2, label='Median')\n", + "#plt.legend(title='%i Filter Changes'%(len(index_filter_change_close_loop)))\n", + "plt.legend(title='%i Filter Changes'%(len(x_array)))\n", + "\n", + "plt.xlabel('Visits after filter change (zero indexed)')\n", + "plt.ylabel(metric)\n", + "#plt.ylabel('Change in aos_fwhm relative to visit before filter change')\n", + "#plt.ylim(0., 1.)\n", + "#plt.ylim(-1.2, 1.2)\n", + "plt.ylim(0.5, 1.5)\n", + "if include is None and exclude is None:\n", + " plt.title(f\"{day_obs_min}-{day_obs_max}\")\n", + "elif include is not None:\n", + " plt.title(f\"{day_obs_min}-{day_obs_max}: include to {include} OR from {include}\")\n", + "else:\n", + " plt.title(f\"{day_obs_min}-{day_obs_max}: exclude to {exclude} OR from {exclude}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6e6c0ad2-4ee4-425a-842f-e515b70b5c92", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e8628cdc-4d8c-4f0b-93a3-00f8ec8f898a", + "metadata": {}, + "outputs": [], + "source": [ + "metric = 'aos_fwhm'\n", + "#metric = 'z4_fwhm'\n", + "#metric = 'psf_fwhm_50'\n", + "\n", + "before_filter = visits_summary['band'].iloc[index_filter_change_close_loop - 1].values\n", + "after_filter = visits_summary['band'].iloc[index_filter_change_close_loop].values\n", + "\n", + "selection = (before_filter == 'y') | (after_filter == 'y')\n", + "\n", + "visits_summary[metric].iloc[index_filter_change_close_loop]\n", + "\n", + "bins = np.linspace(0.3, 1.1, 41)\n", + "#bins = np.linspace(-1., 1., 41)\n", + "#bins = np.linspace(0.5, 1.5, 41)\n", + "\n", + "plt.figure()\n", + "plt.hist(visits_summary[metric].iloc[index_filter_change_close_loop][selection], bins=bins, histtype='step', color='red', label='Filter change involves y')\n", + "plt.hist(visits_summary[metric].iloc[index_filter_change_close_loop][~selection], bins=bins, histtype='step', color='black', label='Filter change does NOT involve y')\n", + "plt.legend()\n", + "plt.xlabel(metric)\n", + "plt.ylabel('Number of Visits')\n", + "\n", + "bins = np.linspace(-0.6, 0.6, 21)\n", + "\n", + "plt.figure()\n", + "\n", + "x = visits_summary[metric].iloc[index_filter_change_close_loop].values - visits_summary[metric].iloc[index_filter_change_close_loop - 1].values\n", + "plt.hist(x[selection], \n", + " bins=bins, histtype='step', color='red', label='to y OR from y')\n", + "plt.hist(x[~selection], \n", + " bins=bins, histtype='step', color='black', label='other')\n", + "\n", + "plt.axvline(np.median(x[selection]), c='red', ls='--')\n", + "plt.axvline(np.median(x[~selection]), c='black', ls='--')\n", + "\n", + "plt.legend()\n", + "plt.xlabel(f\"$\\Delta$ {metric}\")\n", + "plt.ylabel('Number of Visits')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "374cfa5a-a108-43cb-95b4-74e06170b5c9", + "metadata": {}, + "outputs": [], + "source": [ + "visits_summary[metric].iloc[index_filter_change_close_loop].values - visits_summary[metric].iloc[index_filter_change_close_loop - 1].values\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "73eedf60-b507-4c7b-8672-adbb443a01bb", + "metadata": {}, + "outputs": [], + "source": [ + "len(before_filter)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f5c6c2d0-6e21-45ca-8ac3-07b7beaa898f", + "metadata": {}, + "outputs": [], + "source": [ + "np.arange(-10, nvisits)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "af1fcd75-db61-47de-8c97-695db86781bd", + "metadata": {}, + "outputs": [], + "source": [ + "((1. / 2068) - (1. / 3397)) / (1. / 2068)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "361fffdb-f62d-44ea-a5a2-977bb1c4d853", + "metadata": {}, + "outputs": [], + "source": [ + "(1. / 3397) / (1. / 2068)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "23e50549-ff9b-4214-8ef7-eb3c2a95c671", + "metadata": {}, + "outputs": [], + "source": [ + "np.arange(6)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6bad7b50-8cd5-4323-86be-1e86438b924d", + "metadata": {}, + "outputs": [], + "source": [ + "np.sum(~selection)" + ] + }, + { + "cell_type": "markdown", + "id": "b333e980-36a4-40e6-8093-c958787dde14", + "metadata": { + "jp-MarkdownHeadingCollapsed": true + }, + "source": [ + "## Rubin First Alerts" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "be6e2dc9-30f8-43e5-913f-52b073008530", + "metadata": {}, + "outputs": [], + "source": [ + "selection_first_alerts = (visits_summary['science_program'] == 'BLOCK-421') \\\n", + " & (visits_summary['target_name'].str.contains('m49') \\\n", + " | visits_summary['target_name'].str.contains('cosmos') \\\n", + " | visits_summary['target_name'].str.contains('edfs_a') \\\n", + " | visits_summary['target_name'].str.contains('edfs_b') )\n", + "#selection = (visits_summary['science_program'] == 'BLOCK-421')\n", + "print(np.sum(selection_first_alerts))\n", + "visits_summary[selection_first_alerts]['target_name'].value_counts()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ea3c6717-2765-4ad2-ba71-0f432fe6e649", + "metadata": {}, + "outputs": [], + "source": [ + "selection_pre_lsst = (visits_summary['day_obs'] >= 20251230) & (visits_summary['day_obs'] <= 20260127) & (visits_summary['band'] != 'y')\n", + "print(np.sum(selection_pre_lsst))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "05b4072d-2b21-4a75-93a9-9161b9a0bc01", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "\n", + "bins = np.linspace(0.2, 1.0, 40)\n", + "\n", + "plt.hist(visits_summary[selection_first_alerts]['aos_fwhm'], histtype='step', bins=bins, density=True, cumulative=True, lw=2, \n", + " label='First Alerts Targets\\n(EDFS, COSMOS, M49)\\n(%i visits)'%(np.sum(selection_first_alerts)))\n", + "plt.hist(visits_summary[selection_pre_lsst]['aos_fwhm'], histtype='step', bins=bins, density=True, cumulative=True, lw=2, \n", + " label='Standard Pre-LSST\\n(%i visits)'%(np.sum(selection_pre_lsst)))\n", + "plt.legend()\n", + "\n", + "plt.xlabel('AOS FWHM (arcsec)')\n", + "plt.ylabel('CDF')\n", + "plt.xlim(0.2, 1.0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6e4baa2d-00f5-4cd9-be2d-cf417fe418fb", + "metadata": {}, + "outputs": [], + "source": [ + "# Identify filter changes\n", + "selection_filter_change = (visits_summary[\"band\"][1:].values != visits_summary[\"band\"][0:-1].values)\n", + "selection_consecutive = np.diff(visits_summary['seq_num']) == 1\n", + "selection = (selection_consecutive & selection_filter_change & selection_first_alerts[1:])\n", + "\n", + "index_filter_change = np.nonzero(selection)[0] + 1\n", + "\n", + "index_filter_change" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8ca4e2f6-81bb-412d-b283-92cacd2b6962", + "metadata": {}, + "outputs": [], + "source": [ + "for index in index_filter_change:\n", + " print(visits_summary['band'].values[index - 1], \n", + " visits_summary['band'].values[index], \n", + " visits_summary['science_program'].values[index], \n", + " visits_summary['target_name'].values[index])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2c1ec070-3474-4b17-8f2a-5b3df8d97413", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "\n", + "nvisits = 15\n", + "\n", + "metric = 'aos_fwhm'\n", + "#metric = 'z4_fwhm'\n", + "#metric = 'psf_fwhm_50'\n", + "#metric = 'psf_fwhm_zenith_500nm_50'\n", + "\n", + "x_array = []\n", + "y_array = []\n", + "\n", + "for ii in range(0, len(index_filter_change)):\n", + " x = np.arange(nvisits) + np.random.uniform(-0.1, 0.1)\n", + " y = visits_summary[metric].iloc[index_filter_change[ii] + np.arange(nvisits)]\n", + " plt.plot(x, y, lw=1)\n", + "\n", + " x_array.append(x)\n", + " y_array.append(y.values)\n", + "\n", + "\n", + "x_median = []\n", + "y_median = []\n", + "for ii in range(0, nvisits):\n", + " y = []\n", + " for jj in range(0, len(y_array)):\n", + " if len(y_array[jj]) >= ii:\n", + " y.append(y_array[jj][ii])\n", + " \n", + " if len(x) >= 3:\n", + " x_median.append(ii)\n", + " y_median.append(np.nanmedian(y))\n", + " \n", + "plt.plot(x_median, y_median, c='black', lw=2, label='Median')\n", + "plt.xlabel('Visits after filter change (zero indexed)')\n", + "plt.ylabel(metric)\n", + "#plt.ylim(0.4, 0.5)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d3a88bf5-7e9e-44e1-a730-db30e91cfec2", + "metadata": {}, + "outputs": [], + "source": [ + "values = []\n", + "for y in y_array:\n", + " values.append(y[6:])\n", + "values = np.concatenate(values)\n", + "\n", + "plt.figure()\n", + "plt.hist(values)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "81180723-29fb-4810-a295-eb6e96581f31", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "\n", + "bins = np.linspace(0.2, 1.0, 40)\n", + "\n", + "cumulative=True\n", + "\n", + "plt.hist(visits_summary[selection_first_alerts]['aos_fwhm'], histtype='stepfilled', bins=bins, density=True, cumulative=cumulative, lw=1, color='red', alpha=0.2,\n", + " label='First Alerts Targets\\n(EDFS, COSMOS, M49)\\n(%i visits)'%(np.sum(selection_first_alerts)))\n", + "plt.hist(values, histtype='step', bins=bins, density=True, cumulative=cumulative, lw=2, color='red', ls='-',\n", + " label='First Alerts Targets\\n(EDFS, COSMOS, M49)\\nAt least 6 visits after filter change')\n", + "plt.hist(visits_summary[selection_pre_lsst]['aos_fwhm'], histtype='step', bins=bins, density=True, cumulative=cumulative, lw=2, color='black',\n", + " label='Standard Pre-LSST\\n(%i visits)'%(np.sum(selection_pre_lsst)))\n", + "plt.legend()\n", + "\n", + "plt.xlabel('AOS FWHM (arcsec)')\n", + "plt.ylabel('CDF')\n", + "plt.xlim(0.2, 1.0)\n", + "\n", + "print(np.nanmedian(visits_summary[selection_first_alerts]['aos_fwhm']))\n", + "print(np.median(values))\n", + "print(np.nanmedian(visits_summary[selection_pre_lsst]['aos_fwhm']))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b9c037c5-8851-4d15-8fac-ef7f653892c4", + "metadata": {}, + "outputs": [], + "source": [ + "from astropy.time import Time\n", + "\n", + "first_alert_day = Time(61094.5, format='mjd')\n", + "first_alert_day.iso" + ] + }, + { + "cell_type": "markdown", + "id": "dcfb8059-98d0-402a-9e6c-01750cbf854c", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-09T21:37:05.488272Z", + "iopub.status.busy": "2026-03-09T21:37:05.488044Z", + "iopub.status.idle": "2026-03-09T21:37:05.490791Z", + "shell.execute_reply": "2026-03-09T21:37:05.490170Z", + "shell.execute_reply.started": "2026-03-09T21:37:05.488255Z" + }, + "jp-MarkdownHeadingCollapsed": true + }, + "source": [ + "# PSF Characterization" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "868991a7-249e-4f80-92e1-0a733ee970b6", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "plt.hexbin(ccdvisits[\"psf_fwhm_area\"], ccdvisits[\"psf_fwhm\"])\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2e10ad0c-9a0f-4c0d-94b7-8d34c62ffd6b", + "metadata": {}, + "outputs": [], + "source": [ + "#value = ccdvisits[\"psf_fwhm_area\"] / ccdvisits[\"psf_fwhm\"]\n", + "value = visits_summary['psf_fwhm_area_50'] / visits_summary['psf_fwhm_50']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cf2d593b-e364-4212-8737-6249f6332093", + "metadata": {}, + "outputs": [], + "source": [ + "#selection = (ccdvisits[\"donut_blur_fwhm\"] < 0.7) \n", + "selection = (visits_summary[\"donut_blur_fwhm\"] < 0.8) \n", + "#selection = (ccdvisits[\"band\"] != 'y')\n", + "#selection = (visits_summary[\"band\"] != 'y')\n", + "\n", + "plt.figure()\n", + "#plt.hist(value, bins=np.linspace(0.9, 1.5, 101), density=True, histtype='step', label='Selection')\n", + "plt.hist(value[selection], bins=np.linspace(0.9, 1.5, 101), density=True, histtype='step', label='Selection')\n", + "plt.hist(value[~selection], bins=np.linspace(0.9, 1.5, 101), density=True, histtype='step', label='NOT Selection')\n", + "plt.yscale('log')\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "80329d80-fe28-4f21-bfb3-e4d871fd9728", + "metadata": {}, + "outputs": [], + "source": [ + "selection = ((visits_summary['psf_fwhm_area_50'] / visits_summary['psf_fwhm_50']) > 1.25) & (visits_summary['day_obs'] > 20260301)\n", + "\n", + "visits_summary[['day_obs', 'seq_num']][selection]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3dc0f1c2-148e-4cde-97df-ffb8f33f1682", + "metadata": {}, + "outputs": [], + "source": [ + "selection = (visits_summary['day_obs'] == 20260303) & (visits_summary['seq_num'] == 343)\n", + "visits_summary[['psf_fwhm_area_50', 'psf_fwhm_50']][selection]" + ] + }, + { + "cell_type": "markdown", + "id": "9499a25c-8f89-4ccb-a06a-c3645fdfa25c", + "metadata": { + "jp-MarkdownHeadingCollapsed": true + }, + "source": [ + "# Individual Day Summary" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "518b35db-63f9-40d2-b14a-14cfc5187aef", + "metadata": {}, + "outputs": [], + "source": [ + "day_obs = 20260414\n", + "selection = (visits_summary['day_obs'] == day_obs)\n", + "\n", + "print(np.sum(selection))\n", + "\n", + "fig, ax = plt.subplots(3, sharex=True, figsize=(16, 10))\n", + "\n", + "ax[0].errorbar(visits_summary['seq_num'][selection], \n", + " visits_summary['psf_fwhm_50'][selection],\n", + " yerr=np.array([\n", + " visits_summary['psf_fwhm_50'][selection] - visits_summary['psf_fwhm_05'][selection], \n", + " visits_summary['psf_fwhm_95'][selection] - visits_summary['psf_fwhm_50'][selection],\n", + " ]), \n", + " fmt='_', label='PSF FWHM: 5, 50, 95 percentile')\n", + "ax[0].errorbar(visits_summary['seq_num'][selection], \n", + " visits_summary['psf_e_50'][selection],\n", + " yerr=np.array([\n", + " visits_summary['psf_e_50'][selection] - visits_summary['psf_e_05'][selection], \n", + " visits_summary['psf_e_95'][selection] - visits_summary['psf_e_50'][selection],\n", + " ]), \n", + " fmt='_', label='PSF e: 5, 50, 95 percentile')\n", + "#ax[0].errorbar(visits_summary['seq_num'][selection] + 0.2, \n", + "# 0.2 * visits_summary['psf_e_unnormalized_50'][selection],\n", + "# yerr=np.array([\n", + "# 0.2 * (visits_summary['psf_e_unnormalized_50'][selection] - visits_summary['psf_e_unnormalized_05'][selection]), \n", + "# 0.2 * (visits_summary['psf_e_unnormalized_95'][selection] - visits_summary['psf_e_unnormalized_50'][selection]),\n", + "# ]), \n", + "# fmt='_', label='PSF e (unnormalized): 5, 50, 95 percentile', c='magenta')\n", + "#ax[0].scatter(visits_summary['seq_num'], visits_summary['psf_fwhm_95'])\n", + "#ax[0].scatter(visits_summary['seq_num'], visits_summary['psf_fwhm_50'])\n", + "#ax[0].scatter(visits_summary['seq_num'], visits_summary['psf_fwhm_05'])\n", + "ax[0].scatter(visits_summary['seq_num'][selection], visits_summary['psf_fwhm_95_05'][selection], c='red', marker='x', label='sqrt(fwhm_95^2 - fwhm_5^2)')\n", + "ax[0].scatter(visits_summary['seq_num'][selection], visits_summary['donut_blur_fwhm'][selection], c='black', marker='o', label='Donut Blur')\n", + "ax[0].scatter(visits_summary['seq_num'][selection], visits_summary['psf_fwhm_model'][selection], c='0.8', marker='o', label='Model')\n", + "ax[0].scatter(visits_summary['seq_num'][selection], visits_summary['aos_fwhm'][selection], c='green', marker='+', label='AOS Residual')\n", + "#ax[0].scatter(visits_summary['seq_num'][selection], np.fabs(visits_summary[\"z4_fwhm\"][selection] + 0.25), c='magenta', marker='o', label='z4_fwhm')\n", + "ax[0].set_ylabel('Delivered Image Quality Metric\\n(arcsec)')\n", + "#ax[0].scatter(visit_seq_num, np.sqrt(visit_donut_blur_fwhm**2 + visit_psf_fwhm_95_05**2), c='green', marker='+', label='AOS Residual')\n", + "ax[0].legend(loc='upper right')\n", + "#ax[0].set_ylim(0., 2.0)\n", + "ax[0].set_ylim(0., 1.4)\n", + "ax[0].grid(True)\n", + "\n", + "for band in bands:\n", + " band_selection = selection & (visits_summary[\"band\"] == band)\n", + " ax[1].scatter(visits_summary['seq_num'][band_selection], visits_summary['altitude'][band_selection], marker='.', s=50, c=colors[band], label=band)\n", + "#ax[1].scatter(visits_summary['seq_num'], visits_summary['altitude'], marker='.', s=10, c='black')\n", + "ax[1].set_ylabel('Elevation\\n(deg)')\n", + "ax[1].grid(True)\n", + "ax[1].legend()\n", + "\n", + "for band in bands:\n", + " band_selection = selection & (visits_summary[\"band\"] == band)\n", + " ax[2].scatter(visits_summary['seq_num'][band_selection], visits_summary['physical_rotator_angle'][band_selection], marker='.', s=50, c=colors[band], label=band)\n", + "#ax[2].scatter(visits_summary['seq_num'], visits_summary['physical_rotator_angle'], marker='.', s=10, c='black')\n", + "ax[2].set_ylabel('Phys Rot Angle\\n(deg)')\n", + "ax[2].grid(True)\n", + "\n", + "ax[-1].set_xlabel('seq_num')\n", + "\n", + "plt.suptitle('day_obs = %s'%(day_obs))\n", + "\n", + "#plt.suptitle('day_obs = %s'%(day_obs_min))\n", + "plt.subplots_adjust(hspace=0)\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a290968d-6e25-4dea-95f6-7db278077564", + "metadata": {}, + "outputs": [], + "source": [ + "selection = np.isin(visits_summary['day_obs'], [20251026])\n", + "\n", + "bins = np.linspace(0., 1.5, 40)\n", + "\n", + "plt.figure()\n", + "#plt.hist(visits_summary['aos_fwhm'], bins=bins, density=True, histtype='step')\n", + "#plt.hist(visits_summary['aos_fwhm'][selection], bins=bins, density=True, histtype='step')\n", + "plt.hist(visits_summary['psf_fwhm_95_05'], bins=bins, density=True, histtype='step', lw=2, label='All (nvisits = %i)'%(len(selection)))\n", + "plt.hist(visits_summary['psf_fwhm_95_05'][selection], bins=bins, density=True, histtype='step', lw=2, label='20251026 (nvisits = %i)'%(np.sum(selection)))\n", + "plt.xlabel('sqrt(fwhm_95^2 - fwhm_5^2) (arcsec)')\n", + "plt.xlim(0., 1.5)\n", + "plt.legend()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a53a91ec-6956-4598-a713-c561c406cd41", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "plt.scatter(visits_summary['aos_fwhm'], visits_summary['psf_fwhm_95_05'], s=1)\n", + "plt.xlim(0.2, 1.2)\n", + "plt.ylim(0.2, 1.2)" + ] + }, + { + "cell_type": "markdown", + "id": "6dd8e00e-609e-42b4-b012-ea5933b61b51", + "metadata": {}, + "source": [ + "# Day Obs Summary" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9dc096cd-08ef-45c5-b3da-656a498cdca1", + "metadata": {}, + "outputs": [], + "source": [ + "#for key in visits_summary.keys():\n", + "# visits_summary[key] = pd.to_numeric(visits_summary[key], errors=\"coerce\")\n", + "\n", + "# Exclude y due to issues with estimating donut blur\n", + "selection = (visits_summary['band'] != \"y\") & ~visits_summary['observation_reason'].str.contains('filter_change_close_loop')\n", + "\n", + "groups = visits_summary[selection].groupby('day_obs')\n", + "day_obs_summary = pd.DataFrame({\n", + " 'day_obs': groups['day_obs'].first(),\n", + " 'n_visits': groups['day_obs'].count(),\n", + " 'psf_fwhm_95_05_low': groups['psf_fwhm_95_05'].quantile(0.10),\n", + " 'psf_fwhm_95_05_50': groups['psf_fwhm_95_05'].quantile(0.50),\n", + " 'psf_fwhm_95_05_high': groups['psf_fwhm_95_05'].quantile(0.90),\n", + " 'aos_fwhm_low': groups['aos_fwhm'].quantile(0.10),\n", + " 'aos_fwhm_50': groups['aos_fwhm'].quantile(0.50),\n", + " 'aos_fwhm_high': groups['aos_fwhm'].quantile(0.90),\n", + " 'aos_cam_fwhm_low': groups['aos_cam_fwhm'].quantile(0.10),\n", + " 'aos_cam_fwhm_50': groups['aos_cam_fwhm'].quantile(0.50),\n", + " 'aos_cam_fwhm_high': groups['aos_cam_fwhm'].quantile(0.90),\n", + " 'sys_50_low': groups['sys_50'].quantile(0.10),\n", + " 'sys_50_50': groups['sys_50'].quantile(0.50),\n", + " 'sys_50_high': groups['sys_50'].quantile(0.90),\n", + " 'psf_e_50_low': groups['psf_e_50'].quantile(0.10),\n", + " 'psf_e_50_50': groups['psf_e_50'].quantile(0.50),\n", + " 'psf_e_50_high': groups['psf_e_50'].quantile(0.90),\n", + " 'psf_e1_50_low': groups['psf_e1_50'].quantile(0.10),\n", + " 'psf_e1_50_50': groups['psf_e1_50'].quantile(0.50),\n", + " 'psf_e1_50_high': groups['psf_e1_50'].quantile(0.90),\n", + " 'psf_e2_50_low': groups['psf_e2_50'].quantile(0.10),\n", + " 'psf_e2_50_50': groups['psf_e2_50'].quantile(0.50),\n", + " 'psf_e2_50_high': groups['psf_e2_50'].quantile(0.90),\n", + " 'psf_fwhm_50_low': groups['psf_fwhm_50'].quantile(0.10),\n", + " 'psf_fwhm_50_50': groups['psf_fwhm_50'].quantile(0.50),\n", + " 'psf_fwhm_50_high': groups['psf_fwhm_50'].quantile(0.90),\n", + " 'donut_blur_atm_fwhm_low': groups['donut_blur_atm_fwhm'].quantile(0.10),\n", + " 'donut_blur_atm_fwhm_50': groups['donut_blur_atm_fwhm'].quantile(0.50),\n", + " 'donut_blur_atm_fwhm_high': groups['donut_blur_atm_fwhm'].quantile(0.90),\n", + "})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7d9cfbeb-2674-4497-b84b-0e0fc67dd94f", + "metadata": {}, + "outputs": [], + "source": [ + "day_obs_summary " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "628453fb-b3b8-4563-89d4-bd99d903e1f5", + "metadata": {}, + "outputs": [], + "source": [ + "#day_obs_summary['psf_fwhm_50_50']\n", + "#day_obs_summary['aos_cam_fwhm_50']\n", + "day_obs_summary['psf_e_50_50']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8674f49e-8385-47f3-8aca-6d1e02b0e783", + "metadata": {}, + "outputs": [], + "source": [ + "selection = np.tile(True, len(day_obs_summary))\n", + "\n", + "xticks = np.arange(len(day_obs_summary[selection]))\n", + "\n", + "plt.figure(figsize=(10, 6))\n", + "plt.errorbar(xticks - 0.15, day_obs_summary['psf_fwhm_95_05_50'][selection],\n", + " yerr=np.array([day_obs_summary['psf_fwhm_95_05_50'][selection] - day_obs_summary['psf_fwhm_95_05_low'][selection], \n", + " day_obs_summary['psf_fwhm_95_05_high'][selection] - day_obs_summary['psf_fwhm_95_05_50'][selection]]), \n", + " fmt='_', c='tab:blue', label='sqrt(fwhm_95^2 - fwhm_5^2): 10th, 50th, 90th Percentile')\n", + "#plt.errorbar(xticks + 0.2, day_obs_summary['aos_fwhm_50'][selection],\n", + "# yerr=np.array([day_obs_summary['aos_fwhm_50'][selection] - day_obs_summary['aos_fwhm_low'][selection], \n", + "# day_obs_summary['aos_fwhm_high'][selection] - day_obs_summary['aos_fwhm_50'][selection]]), \n", + "# fmt='_', c='magenta', label='aos_fwhm: 10th, 50th, 90th Percentile')\n", + "plt.errorbar(xticks + 0., day_obs_summary['aos_cam_fwhm_50'][selection],\n", + " yerr=np.array([day_obs_summary['aos_cam_fwhm_50'][selection] - day_obs_summary['aos_cam_fwhm_low'][selection], \n", + " day_obs_summary['aos_cam_fwhm_high'][selection] - day_obs_summary['aos_cam_fwhm_50'][selection]]), \n", + " fmt='_', c='tab:purple', label='sqrt(aos_fwhm^2 + cam_fwhm^2): 10th, 50th, 90th Percentile')\n", + "plt.errorbar(xticks + 0.15, day_obs_summary['sys_50_50'][selection],\n", + " yerr=np.array([day_obs_summary['sys_50_50'][selection] - day_obs_summary['sys_50_low'][selection], \n", + " day_obs_summary['sys_50_high'][selection] - day_obs_summary['sys_50_50'][selection]]), \n", + " fmt='_', c='tab:red', label='sqrt(fwhm_50^2 - (donut_blur^2 - cam_fwhm^2)): 10th, 50th, 90th Percentile')\n", + "plt.errorbar(xticks + 0., day_obs_summary['psf_e_50_50'][selection],\n", + " yerr=np.array([day_obs_summary['psf_e_50_50'][selection] - day_obs_summary['psf_e_50_low'][selection], \n", + " day_obs_summary['psf_e_50_high'][selection] - day_obs_summary['psf_e_50_50'][selection]]), \n", + " fmt='_', c='tab:orange', label='psf_e: 10th, 50th, 90th Percentile')\n", + "\n", + "\n", + "plt.errorbar(xticks - 0.075, day_obs_summary['psf_fwhm_50_50'][selection],\n", + " yerr=np.array([day_obs_summary['psf_fwhm_50_50'][selection] - day_obs_summary['psf_fwhm_50_low'][selection], \n", + " day_obs_summary['psf_fwhm_50_high'][selection] - day_obs_summary['psf_fwhm_50_50'][selection]]), \n", + " fmt='_', c='black', label='fwhm_50: 10th, 50th, 90th Percentile')\n", + "plt.errorbar(xticks + 0.075, day_obs_summary['donut_blur_atm_fwhm_50'][selection],\n", + " yerr=np.array([day_obs_summary['donut_blur_atm_fwhm_50'][selection] - day_obs_summary['donut_blur_atm_fwhm_low'][selection], \n", + " day_obs_summary['donut_blur_atm_fwhm_high'][selection] - day_obs_summary['donut_blur_atm_fwhm_50'][selection]]), \n", + " fmt='_', c='0.5', label='sqrt(donut_blur^2 - cam_fwhm^2): 10th, 50th, 90th Percentile')\n", + "\n", + "plt.legend()\n", + "plt.xticks(xticks, day_obs_summary['day_obs'][selection], rotation=90)\n", + "plt.ylim(0, plt.ylim()[-1])\n", + "#for y in [0.04, 0.2, 0.4, 0.6]:\n", + "# plt.axhline(y, c='0.75', ls='--', zorder=0)\n", + "\n", + "plt.axhline(1., c='0.25', ls=':', zorder=0)\n", + "plt.axhline(0.45, c='red', ls=':', zorder=0)\n", + "plt.axhline(0.04, c='orange', ls=':', zorder=0)\n", + "\n", + "plt.ylabel('Image Quality Contribution (arcsec)')\n", + "plt.tight_layout()\n", + "#plt.savefig('system_diq_timeline.pdf')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "69ac1f5a-45f2-4e19-b54a-150f735e441e", + "metadata": {}, + "outputs": [], + "source": [ + "from astropy.time import Time\n", + "\n", + "def dayObsToTime(day_obs):\n", + " year = str(day_obs)[0:4]\n", + " month = str(day_obs)[4:6]\n", + " day = str(day_obs)[6:8]\n", + " time = Time('%s-%s-%s 12:00:00'%(year, month, day), format='iso')\n", + " return time\n", + "\n", + "selection = np.tile(True, len(day_obs_summary))\n", + "\n", + "#xticks = np.arange(len(day_obs_summary[selection]))\n", + "xticks = np.array([dayObsToTime(_).mjd for _ in day_obs_summary['day_obs'][selection]])\n", + "\n", + "print(xticks)\n", + "\n", + "fig, ax = plt.subplots(figsize=(15, 6), dpi=200)\n", + "a = np.linspace(0., 1., 10)\n", + "\n", + "variation = ax.errorbar(\n", + " xticks - 0.15, day_obs_summary['psf_fwhm_95_05_50'][selection],\n", + " yerr=np.array([\n", + " day_obs_summary['psf_fwhm_95_05_50'][selection] - day_obs_summary['psf_fwhm_95_05_low'][selection], \n", + " day_obs_summary['psf_fwhm_95_05_high'][selection] - day_obs_summary['psf_fwhm_95_05_50'][selection]\n", + " ]), \n", + " fmt='_', c='tab:blue', label='Variation Across FoV'\n", + " #label='sqrt(fwhm_95^2 - fwhm_5^2)'\n", + ")\n", + "aos = ax.errorbar(\n", + " xticks + 0., day_obs_summary['aos_cam_fwhm_50'][selection],\n", + " yerr=np.array([day_obs_summary['aos_cam_fwhm_50'][selection] - day_obs_summary['aos_cam_fwhm_low'][selection], \n", + " day_obs_summary['aos_cam_fwhm_high'][selection] - day_obs_summary['aos_cam_fwhm_50'][selection]]), \n", + " fmt='_', c='tab:purple', label='Optics + Camera Diffusion'\n", + " #label='sqrt(aos_fwhm^2 + cam_fwhm^2)'\n", + ")\n", + "sys = ax.errorbar(\n", + " xticks + 0.15, day_obs_summary['sys_50_50'][selection],\n", + " yerr=np.array([day_obs_summary['sys_50_50'][selection] - day_obs_summary['sys_50_low'][selection], \n", + " day_obs_summary['sys_50_high'][selection] - day_obs_summary['sys_50_50'][selection]]), \n", + " fmt='_', c='tab:red', label='Delivered - Estimated Atmosphere'\n", + " #label='sqrt(fwhm_50^2 - (donut_blur^2 - cam_fwhm^2))'\n", + ")\n", + "\n", + "ellipticity = ax.errorbar(\n", + " xticks + 0., day_obs_summary['psf_e_50_50'][selection],\n", + " yerr=np.array([day_obs_summary['psf_e_50_50'][selection] - day_obs_summary['psf_e_50_low'][selection], \n", + " day_obs_summary['psf_e_50_high'][selection] - day_obs_summary['psf_e_50_50'][selection]]), \n", + " fmt='_', c='tab:orange', label='PSF Ellipticity'\n", + " #label='psf_e'\n", + ")\n", + "\"\"\"\n", + "ellipticity_1 = ax.errorbar(\n", + " xticks - 0.15, day_obs_summary['psf_e1_50_50'][selection],\n", + " yerr=np.array([day_obs_summary['psf_e1_50_50'][selection] - day_obs_summary['psf_e1_50_low'][selection], \n", + " day_obs_summary['psf_e1_50_high'][selection] - day_obs_summary['psf_e1_50_50'][selection]]), \n", + " fmt='_', c='tab:orange',\n", + " #label='psf_e'\n", + ")\n", + "ellipticity_2 = ax.errorbar(\n", + " xticks + 0.15, day_obs_summary['psf_e2_50_50'][selection],\n", + " yerr=np.array([day_obs_summary['psf_e2_50_50'][selection] - day_obs_summary['psf_e2_50_low'][selection], \n", + " day_obs_summary['psf_e2_50_high'][selection] - day_obs_summary['psf_e2_50_50'][selection]]), \n", + " fmt='_', c='tab:orange',\n", + " #label='psf_e'\n", + ")\n", + "\"\"\"\n", + "\n", + "psf = ax.errorbar(\n", + " xticks - 0.075, day_obs_summary['psf_fwhm_50_50'][selection],\n", + " yerr=np.array([day_obs_summary['psf_fwhm_50_50'][selection] - day_obs_summary['psf_fwhm_50_low'][selection], \n", + " day_obs_summary['psf_fwhm_50_high'][selection] - day_obs_summary['psf_fwhm_50_50'][selection]]), \n", + " fmt='_', c='black', label='Delivered Median PSF FWHM' \n", + " #label='fwhm_50'\n", + ")\n", + "atm = ax.errorbar(\n", + " xticks + 0.075, day_obs_summary['donut_blur_atm_fwhm_50'][selection],\n", + " yerr=np.array([day_obs_summary['donut_blur_atm_fwhm_50'][selection] - day_obs_summary['donut_blur_atm_fwhm_low'][selection], \n", + " day_obs_summary['donut_blur_atm_fwhm_high'][selection] - day_obs_summary['donut_blur_atm_fwhm_50'][selection]]), \n", + " fmt='_', c='0.5', label='Estimated Atmosphere'\n", + " #label='sqrt(donut_blur^2 - cam_fwhm^2)'\n", + ")\n", + "\n", + "\n", + "first_legend = ax.legend(handles=[variation, aos, sys], loc='upper left', title='Estimated System Contribution')\n", + "#first_legend = ax.legend(handles=[variation, aos], loc='upper left', title='Estimated System Contribution')\n", + "second_legend = ax.legend(handles=[psf, atm, ellipticity], loc='upper right', title='Measured Quantities')\n", + "ax.add_artist(first_legend)\n", + "\n", + "ax.set_xticks(xticks, day_obs_summary['day_obs'][selection], rotation=90, fontsize=7.5)\n", + "#ax.set_ylim(0, plt.ylim()[-1])\n", + "#ax.set_ylim(0, 2.5)\n", + "ax.set_ylim(0., 2.0)\n", + "#ax.set_ylim(0.3, 0.5)\n", + "#for y in [0.04, 0.2, 0.4, 0.6]:\n", + "# plt.axhline(y, c='0.75', ls='--', zorder=0)\n", + "\n", + "#ax.axhline(1., c='0.25', ls=':', zorder=0)\n", + "ax.axhline(0.45, c='red', ls=':', zorder=0)\n", + "ax.axhline(0.04, c='orange', ls=':', zorder=0)\n", + "\n", + "ax.set_ylabel('Image Quality Metric (arcsec)')\n", + "plt.title('10th, 50th, 90th Percentiles of Per-visit Quantities for Each Night')\n", + "plt.grid()\n", + "fig.tight_layout()\n", + "#plt.show()\n", + "#plt.savefig('system_diq_timeline.pdf')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "01195e09-9774-4c75-b9e8-425b72261c4b", + "metadata": {}, + "outputs": [], + "source": [ + "xticks = np.arange(len(day_obs_summary))\n", + "\n", + "plt.figure(figsize=(10, 6))\n", + "plt.errorbar(xticks, day_obs_summary['aos_fwhm_50'],\n", + " yerr=np.array([day_obs_summary['aos_fwhm_50'] - day_obs_summary['aos_fwhm_05'], \n", + " day_obs_summary['aos_fwhm_95'] - day_obs_summary['aos_fwhm_50']]), \n", + " fmt='_')\n", + "plt.xticks(xticks, day_obs_summary['day_obs'], rotation=90)\n", + "plt.ylabel('aos_fwhm')\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e8cfa674-7a7b-4954-8bbf-4190a6355558", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "for band in bands:\n", + " selection = (visits_summary[\"band\"] == band)\n", + " plt.scatter(visits_summary['donut_blur_fwhm'][selection], visits_summary['psf_fwhm_50'][selection], s=2, c=colors[band], label=band)\n", + "plt.xlim(0.5, 2.)\n", + "plt.ylim(0.5, 2.)\n", + "plt.legend(markerscale=3)\n", + "plt.xlabel('donut_blur_fwhm (arcsec)')\n", + "plt.ylabel('median_psf_fwhm (arcsec)')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "97fa7ffb-2212-4dd8-a60b-f089f5338b59", + "metadata": {}, + "outputs": [], + "source": [ + "bins=np.linspace(-0.5, 0.5, 50)\n", + "\n", + "plt.figure()\n", + "#for band in bands:\n", + "for band in ['g', 'r', 'i', 'z']:\n", + " selection = (visits_summary[\"band\"] == band) & (visits_summary['psf_fwhm_50'] < 2.0)\n", + " plt.hist(visits_summary['donut_blur_fwhm'][selection] - visits_summary['psf_fwhm_50'][selection], bins=bins, density=True, histtype='step', color=colors[band], label=band, lw=2)\n", + "plt.legend()\n", + "plt.xlabel('median_psf_fwhm - donut_blur_fwhm (arcsec)')\n", + "plt.ylabel('PDF')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4ec2590b-11f5-461b-b16f-c73b4315f8a1", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "plt.scatter(day_obs_summary['donut_blur_atm_fwhm_50'], day_obs_summary['psf_e_50_50'])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "72fdc2f5-c7ad-48b0-a4db-e24d20760eec", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "63c98239-3379-453e-a922-ca285e29355f", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "4605d5f6-2d0a-41ad-8b6d-6ae5a34a326b", + "metadata": { + "jp-MarkdownHeadingCollapsed": true + }, + "source": [ + "# Best Convergence by Day" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ffef8c92-ce8f-44c9-b9e6-b0d81939e0f9", + "metadata": {}, + "outputs": [], + "source": [ + "# Select best visits\n", + "threshold_psf_fwhm_50 = 0.9\n", + "threshold_aos_fwhm = 0.45\n", + "selection = (visits_summary['band'] != \"y\") \\\n", + " & (visits_summary['psf_fwhm_50'] < threshold_psf_fwhm_50) \\\n", + " & (visits_summary['aos_fwhm'] < threshold_aos_fwhm)\n", + "\n", + "groups = visits_summary[selection].groupby('day_obs')\n", + "day_obs_summary = pd.DataFrame({\n", + " 'day_obs': groups['day_obs'].first(),\n", + " 'n_visits': groups['day_obs'].count(),\n", + " 'psf_fwhm_95_05_low': groups['psf_fwhm_95_05'].quantile(0.10),\n", + " 'psf_fwhm_95_05_50': groups['psf_fwhm_95_05'].quantile(0.50),\n", + " 'psf_fwhm_95_05_high': groups['psf_fwhm_95_05'].quantile(0.90),\n", + " 'aos_fwhm_low': groups['aos_fwhm'].quantile(0.10),\n", + " 'aos_fwhm_50': groups['aos_fwhm'].quantile(0.50),\n", + " 'aos_fwhm_high': groups['aos_fwhm'].quantile(0.90),\n", + " 'aos_cam_fwhm_low': groups['aos_cam_fwhm'].quantile(0.10),\n", + " 'aos_cam_fwhm_50': groups['aos_cam_fwhm'].quantile(0.50),\n", + " 'aos_cam_fwhm_high': groups['aos_cam_fwhm'].quantile(0.90),\n", + " 'sys_50_low': groups['sys_50'].quantile(0.10),\n", + " 'sys_50_50': groups['sys_50'].quantile(0.50),\n", + " 'sys_50_high': groups['sys_50'].quantile(0.90),\n", + " 'psf_e_50_low': groups['psf_e_50'].quantile(0.10),\n", + " 'psf_e_50_50': groups['psf_e_50'].quantile(0.50),\n", + " 'psf_e_50_high': groups['psf_e_50'].quantile(0.90),\n", + " 'psf_e1_50_low': groups['psf_e1_50'].quantile(0.10),\n", + " 'psf_e1_50_50': groups['psf_e1_50'].quantile(0.50),\n", + " 'psf_e1_50_high': groups['psf_e1_50'].quantile(0.90),\n", + " 'psf_e2_50_low': groups['psf_e2_50'].quantile(0.10),\n", + " 'psf_e2_50_50': groups['psf_e2_50'].quantile(0.50),\n", + " 'psf_e2_50_high': groups['psf_e2_50'].quantile(0.90),\n", + " 'psf_fwhm_50_low': groups['psf_fwhm_50'].quantile(0.10),\n", + " 'psf_fwhm_50_50': groups['psf_fwhm_50'].quantile(0.50),\n", + " 'psf_fwhm_50_high': groups['psf_fwhm_50'].quantile(0.90),\n", + " 'donut_blur_atm_fwhm_low': groups['donut_blur_atm_fwhm'].quantile(0.10),\n", + " 'donut_blur_atm_fwhm_50': groups['donut_blur_atm_fwhm'].quantile(0.50),\n", + " 'donut_blur_atm_fwhm_high': groups['donut_blur_atm_fwhm'].quantile(0.90),\n", + "})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d6503cc9-3511-4a9e-b8a4-8584d04e1771", + "metadata": {}, + "outputs": [], + "source": [ + "day_obs_summary['n_visits']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6c69705d-5dd7-4fcf-9d5a-4cdb7b6f68de", + "metadata": {}, + "outputs": [], + "source": [ + "from astropy.time import Time\n", + "\n", + "def dayObsToTime(day_obs):\n", + " year = str(day_obs)[0:4]\n", + " month = str(day_obs)[4:6]\n", + " day = str(day_obs)[6:8]\n", + " time = Time('%s-%s-%s 12:00:00'%(year, month, day), format='iso')\n", + " return time\n", + "\n", + "#selection = np.tile(True, len(day_obs_summary))\n", + "selection = (day_obs_summary['n_visits'] > 40)\n", + "\n", + "#xticks = np.arange(len(day_obs_summary[selection]))\n", + "xticks = np.array([dayObsToTime(_).mjd for _ in day_obs_summary['day_obs'][selection]])\n", + "\n", + "fig, ax = plt.subplots(figsize=(10, 6), dpi=200)\n", + "\n", + "plt.scatter(xticks, day_obs_summary['psf_e_50_50'][selection], c='tab:red', label='10th Percentile of Per-visit Median Ellipticity')\n", + "plt.scatter(xticks, day_obs_summary['psf_e_50_low'][selection], c='orange', label='10th Percentile of Per-visit Median Ellipticity')\n", + "plt.scatter(xticks, day_obs_summary['psf_fwhm_95_05_low'][selection], c='tab:blue', label='10th Percentile of Per-visit sqrt(fwhm_95^2 - fwhm_5^2) (arcsec)')\n", + "\n", + "plt.legend()\n", + "\n", + "ax.set_xticks(xticks, day_obs_summary['day_obs'][selection], rotation=90, fontsize=6)\n", + "#ax.set_ylim(0., 0.2)\n", + "ax.axhline(0.04, c='orange', ls=':', zorder=0)\n", + "\n", + "ax.set_ylabel('Image Quality Metric')\n", + "#ax.set_ylabel('10th Percentile of Per-visit Median Ellipticity')\n", + "#ax.set_ylabel('10th Percentile of Per-visit Median Ellipticity')\n", + "#plt.title('10th, 50th, 90th Percentiles of Per-visit Quantities for Each Night')\n", + "plt.title('BLOCK-T539, BLOCK-T698, BLOCK-T641;\\nConsider Only visits with PSF FHWM < %.2f arcsec & AOS Residual < %.2f arcsec'%(threshold_psf_fwhm_50, threshold_aos_fwhm))\n", + "plt.grid()\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "id": "d3ee9406-e927-435a-83e8-687f4e8cf390", + "metadata": { + "jp-MarkdownHeadingCollapsed": true + }, + "source": [ + "# Per sensor" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2492fa24-858c-4646-b7c8-e878e4c9280b", + "metadata": {}, + "outputs": [], + "source": [ + "ccdvisits.columns" + ] + }, + { + "cell_type": "markdown", + "id": "740b9cf4-4f0f-4961-a9d7-40cb93ec61fc", + "metadata": {}, + "source": [ + "- 20260310 (Tue): Scripted AOS closed loop\n", + "- 20260311 (Wed): FBS-driven closed loop w/ MTAOS\n", + "- 20260312 (Thu): FBS-driven closed loop w/ MTAOS: medium gain\n", + "- 20260313 (Fri): FBS-driven closed loop w/ MTAOS: high gain\n", + "- 20260314 (Sat): FBS-driven closed loop w/ MTAOS: DOF\n", + "- 20260318 (Sat): FBS-driven closed loop w/ MTAOS: WCS donut selection\n", + "- 20260325 (): FBS-driven closed loop w/ MTAOS: turn off M2 dz\n", + "- 20260326 (): FBS-driven closed loop w/ MTAOS: TMA fans\n", + "- 20260328 (Sat): FBS-driven closed loop w/ MTAOS: TMA fans\n", + "- 20260402: FBS-driven closed loop w/ MTAOS: simplified control-loop mode discarding the intermediate updates\n", + "- 20260403: FBS-driven closed loop w/ MTAOS: PID coefficient optimization (kp), simplified control-loop mode discarding the intermediate updates\n", + "- 20260403: FBS-driven closed loop w/ MTAOS: KP = 0.3, KI = 0, discard intermediate updates, Danish binning = 2, truncation index = 9 (was 12 on previous nights); intent was to remove spherical vmode" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e7d39e40-e3a8-4fb9-8b14-e112c5a0d6ec", + "metadata": {}, + "outputs": [], + "source": [ + "day_obs_select = 20260414\n", + "threshold_aos_fwhm = 0.45\n", + "selection = (ccdvisits['band'] != 'y') & (ccdvisits['aos_fwhm'] < threshold_aos_fwhm) & (ccdvisits['day_obs'] == day_obs_select) \\\n", + " & (ccdvisits['physical_rotator_angle'] > -10.) & (ccdvisits['physical_rotator_angle'] < 10.)\n", + "\n", + "groups = ccdvisits[selection].groupby('detector')\n", + "detector_summary = pd.DataFrame({\n", + " 'detector': groups['detector'].first(),\n", + " 'n_visits': groups['detector'].count(),\n", + " 'psf_fwhm_05': groups['psf_fwhm'].quantile(0.25),\n", + " 'psf_fwhm_50': groups['psf_fwhm'].quantile(0.50),\n", + " 'psf_fwhm_95': groups['psf_fwhm'].quantile(0.75),\n", + " 'psf_fwhm_area_50': groups['psf_fwhm_area'].quantile(0.50),\n", + " 'psf_e_05': groups['psf_e'].quantile(0.25),\n", + " 'psf_e_50': groups['psf_e'].quantile(0.50),\n", + " 'psf_e_95': groups['psf_e'].quantile(0.75),\n", + " 'psf_e1_05': groups['psf_e1'].quantile(0.25),\n", + " 'psf_e1_50': groups['psf_e1'].quantile(0.50),\n", + " 'psf_e1_95': groups['psf_e1'].quantile(0.75),\n", + " 'psf_e2_05': groups['psf_e2'].quantile(0.25),\n", + " 'psf_e2_50': groups['psf_e2'].quantile(0.50),\n", + " 'psf_e2_95': groups['psf_e2'].quantile(0.75),\n", + "})\n", + "\n", + "\"\"\"\n", + "visits_summary[selection].groupby('day_obs')\n", + "\n", + "\n", + "threshold_psf_fwhm_50 = 0.9\n", + "threshold_aos_fwhm = 0.45\n", + "selection = (visits_summary['band'] != \"y\") \\\n", + " & (visits_summary['psf_fwhm_50'] < threshold_psf_fwhm_50) \\\n", + " & (visits_summary['aos_fwhm'] < threshold_aos_fwhm)\n", + "\n", + "groups = visits_summary[selection].groupby('day_obs')\n", + "day_obs_summary = pd.DataFrame({\n", + " 'day_obs': groups['day_obs'].first(),\n", + " 'n_visits': groups['day_obs'].count(),\n", + " 'psf_fwhm_95_05_low': groups['psf_fwhm_95_05'].quantile(0.10),\n", + " 'psf_fwhm_95_05_50': groups['psf_fwhm_95_05'].quantile(0.50),\n", + " 'psf_fwhm_95_05_high': groups['psf_fwhm_95_05'].quantile(0.90),\n", + " 'aos_fwhm_low': groups['aos_fwhm'].quantile(0.10),\n", + " 'aos_fwhm_50': groups['aos_fwhm'].quantile(0.50),\n", + " 'aos_fwhm_high': groups['aos_fwhm'].quantile(0.90),\n", + " 'aos_cam_fwhm_low': groups['aos_cam_fwhm'].quantile(0.10),\n", + " 'aos_cam_fwhm_50': groups['aos_cam_fwhm'].quantile(0.50),\n", + " 'aos_cam_fwhm_high': groups['aos_cam_fwhm'].quantile(0.90),\n", + " 'sys_50_low': groups['sys_50'].quantile(0.10),\n", + " 'sys_50_50': groups['sys_50'].quantile(0.50),\n", + " 'sys_50_high': groups['sys_50'].quantile(0.90),\n", + " 'psf_e_50_low': groups['psf_e_50'].quantile(0.10),\n", + " 'psf_e_50_50': groups['psf_e_50'].quantile(0.50),\n", + " 'psf_e_50_high': groups['psf_e_50'].quantile(0.90),\n", + " 'psf_e1_50_low': groups['psf_e1_50'].quantile(0.10),\n", + " 'psf_e1_50_50': groups['psf_e1_50'].quantile(0.50),\n", + " 'psf_e1_50_high': groups['psf_e1_50'].quantile(0.90),\n", + " 'psf_e2_50_low': groups['psf_e2_50'].quantile(0.10),\n", + " 'psf_e2_50_50': groups['psf_e2_50'].quantile(0.50),\n", + " 'psf_e2_50_high': groups['psf_e2_50'].quantile(0.90),\n", + " 'psf_fwhm_50_low': groups['psf_fwhm_50'].quantile(0.10),\n", + " 'psf_fwhm_50_50': groups['psf_fwhm_50'].quantile(0.50),\n", + " 'psf_fwhm_50_high': groups['psf_fwhm_50'].quantile(0.90),\n", + " 'donut_blur_atm_fwhm_low': groups['donut_blur_atm_fwhm'].quantile(0.10),\n", + " 'donut_blur_atm_fwhm_50': groups['donut_blur_atm_fwhm'].quantile(0.50),\n", + " 'donut_blur_atm_fwhm_high': groups['donut_blur_atm_fwhm'].quantile(0.90),\n", + "})\n", + "\n", + "groups = ccdvisits.groupby('visit_id')\n", + "visits_summary = pd.DataFrame({\n", + " 'day_obs': groups['day_obs'].first(),\n", + " 'target_name': groups['target_name'].first(),\n", + " 'science_program': groups['science_program'].first(),\n", + " 'observation_reason': groups['observation_reason'].first(),\n", + " 'seq_num': groups['seq_num'].median(),\n", + " 'exp_midpt_mjd': groups['exp_midpt_mjd'].median(),\n", + " 'donut_blur_fwhm': groups['donut_blur_fwhm'].median(),\n", + " 'aos_fwhm': groups['aos_fwhm'].median(),\n", + " 'donut_blur_atm_fwhm': groups['donut_blur_atm_fwhm'].median(),\n", + " 'aos_cam_fwhm': groups['aos_cam_fwhm'].median(),\n", + " 'physical_rotator_angle': groups['physical_rotator_angle'].median(),\n", + " 'altitude': groups['altitude'].median(),\n", + " 'airmass': groups['airmass'].median(),\n", + " 'psf_fwhm_05': groups['psf_fwhm'].quantile(0.05),\n", + " 'psf_fwhm_50': groups['psf_fwhm'].quantile(0.50),\n", + " 'psf_fwhm_95': groups['psf_fwhm'].quantile(0.95),\n", + " 'psf_fwhm_zenith_500nm_50': groups['psf_fwhm_zenith_500nm'].quantile(0.50),\n", + " 'psf_fwhm_area_50': groups['psf_fwhm_area'].quantile(0.50),\n", + " 'psf_e_05': groups['psf_e'].quantile(0.05),\n", + " 'psf_e_50': groups['psf_e'].quantile(0.50),\n", + " 'psf_e_95': groups['psf_e'].quantile(0.95),\n", + " 'psf_e1_05': groups['psf_e1'].quantile(0.05),\n", + " 'psf_e1_50': groups['psf_e1'].quantile(0.50),\n", + " 'psf_e1_95': groups['psf_e1'].quantile(0.95),\n", + " 'psf_e2_05': groups['psf_e2'].quantile(0.05),\n", + " 'psf_e2_50': groups['psf_e2'].quantile(0.50),\n", + " 'psf_e2_95': groups['psf_e2'].quantile(0.95),\n", + " 'psf_e_unnormalized_05': groups['psf_e_unnormalized'].quantile(0.05),\n", + " 'psf_e_unnormalized_50': groups['psf_e_unnormalized'].quantile(0.50),\n", + " 'psf_e_unnormalized_95': groups['psf_e_unnormalized'].quantile(0.95),\n", + " 'band': groups['band'].first(),\n", + "})\n", + "\"\"\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "24b8df85-e0ab-4d30-b425-188590905b9a", + "metadata": {}, + "outputs": [], + "source": [ + "detector_summary" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "afe7ff65-60fe-4355-ad22-f2b4107d5182", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "#plt.scatter(detector_summary['detector'], detector_summary['psf_fwhm_50'])\n", + "\n", + "plt.scatter(detector_summary['detector'], detector_summary['psf_e_50'])\n", + "\n", + "\n", + "\"\"\"\n", + "plt.errorbar(detector_summary['detector'], \n", + " detector_summary['psf_fwhm_50'],\n", + " yerr=np.array([\n", + " detector_summary['psf_fwhm_50'] - detector_summary['psf_fwhm_05'], \n", + " detector_summary['psf_fwhm_95'] - detector_summary['psf_fwhm_50'],\n", + " ]), \n", + " fmt='_', label='PSF FWHM: 5, 50, 95 percentile')\n", + "\"\"\"\n", + "\n", + "\n", + "plt.errorbar(detector_summary['detector'], \n", + " detector_summary['psf_e_50'],\n", + " yerr=np.array([\n", + " detector_summary['psf_e_50'] - detector_summary['psf_e_05'], \n", + " detector_summary['psf_e_95'] - detector_summary['psf_e_50'],\n", + " ]), \n", + " fmt='_', label='PSF e: 5, 50, 95 percentile')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "adc30453-8aaf-49fd-b1a3-b339db87e782", + "metadata": {}, + "outputs": [], + "source": [ + "import lsst.daf.butler\n", + "#import lsst.afw.geom as afwGeom\n", + "import lsst.afw.cameraGeom as cameraGeom\n", + "\n", + "def getCamera(instrument='LSSTCam'):\n", + " repo = '/repo/embargo'\n", + " collection = 'LSSTCam/calib'\n", + " butler = lsst.daf.butler.Butler(repo, collections=collection)\n", + " camera = butler.get('camera', instrument=instrument)\n", + " return camera\n", + "\n", + "def getCorners(det):\n", + " x = np.array([\n", + " det.getBBox().beginX,\n", + " det.getBBox().beginX,\n", + " det.getBBox().endX,\n", + " det.getBBox().endX\n", + " ])\n", + " y = np.array([\n", + " det.getBBox().beginY,\n", + " det.getBBox().endY,\n", + " det.getBBox().endY,\n", + " det.getBBox().beginY,\n", + " ])\n", + " return x, y" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "eb4500e4-98e2-4341-8e5a-65b8563e1800", + "metadata": {}, + "outputs": [], + "source": [ + "camera = getCamera()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e0f21e05-e8f0-43a9-b1e5-773367e33da2", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib as mpl\n", + "import matplotlib.cm as cm" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b134a1c3-0e2d-4135-90dd-f6fc328c95ac", + "metadata": {}, + "outputs": [], + "source": [ + "cmap = cm.get_cmap('viridis')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "173b81e6-4c33-488b-8e47-274d9cb9a1d7", + "metadata": {}, + "outputs": [], + "source": [ + "def rescale(x, minimum, maximum):\n", + " return (x - minimum) / (maximum - minimum)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5bc33e53-8b41-491c-a470-7f186ba0ddc1", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(figsize=(8,6))\n", + "\n", + "#metric = 'psf_fwhm_50'\n", + "metric = 'psf_e_50'\n", + "\n", + "#fc='blue'\n", + "alpha=1\n", + "\n", + "plt.scatter(1.e6 + detector_summary['detector'], 1.e6 + detector_summary['detector'], c=detector_summary[metric])\n", + "\n", + "for det in camera:\n", + "\n", + " if det.getId() in detector_summary['detector']:\n", + "\n", + " value = detector_summary.loc[detector_summary['detector'] == det.getId()][metric].values[0]\n", + " fc = cmap(rescale(\n", + " value, \n", + " np.min(detector_summary[metric]),\n", + " np.max(detector_summary[metric])\n", + " ))\n", + " \n", + " x, y = zip(*[(corner.x, corner.y) for corner in det.getCorners(cameraGeom.FOCAL_PLANE)])\n", + " plt.fill(x, y, fc=fc, alpha=alpha)\n", + "\n", + "plt.colorbar(cmap=cmap.name, label=metric)\n", + "\n", + "plt.xlim(-350, 350)\n", + "plt.ylim(-350, 350)\n", + "\n", + "plt.title('day_obs = %i\\n(n = %i)'%(day_obs_select, np.median(detector_summary['n_visits'])))\n", + "\n", + "\"\"\"\n", + "norm = mpl.colors.Normalize(vmin=0, vmax=1)\n", + "s_mappable = mpl.cm.ScalarMappable(cmap=cmap, norm=norm)\n", + "s_mappable.set_array([])\n", + "cbar = plt.colorbar(s_mappable, cax=plt.gca(), orientation='vertical')\n", + "cbar.set_label('Some relevant units')\n", + "\"\"\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "89a081b3-cb47-443e-917c-73903221f836", + "metadata": {}, + "outputs": [], + "source": [ + "dir(det)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d7f9f604-9d5d-4324-8f62-2eb21c06eadb", + "metadata": {}, + "outputs": [], + "source": [ + "det.getId()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "324ece64-b590-4b6a-bf7f-4a5d3a2008e1", + "metadata": {}, + "outputs": [], + "source": [ + "14 in detector_summary['detector']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c26df2a8-e752-4f2a-af10-5ec0cbede76c", + "metadata": {}, + "outputs": [], + "source": [ + "detector_summary.loc[detector_summary['detector'] == 14]['psf_fwhm_50'].values[0]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c74b80e2-3035-4915-aa84-a8e0f18e783a", + "metadata": {}, + "outputs": [], + "source": [ + "rescale(0.5, 0., 0.5)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0edc4f21-aba9-4b24-9a57-ea9252aa3bbb", + "metadata": {}, + "outputs": [], + "source": [ + "cmap.name" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2f6211c1-1f08-44b5-b75a-c7f9693a9ef3", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(2, 2, figsize=(10, 8), dpi=150)\n", + "\n", + "alpha=1\n", + "\n", + "metric = 'psf_fwhm_50'\n", + "\n", + "scat = ax[0][0].scatter(1.e6 + detector_summary['detector'], 1.e6 + detector_summary['detector'], c=detector_summary[metric])\n", + "for det in camera:\n", + " if det.getId() in detector_summary['detector']:\n", + " value = detector_summary.loc[detector_summary['detector'] == det.getId()][metric].values[0]\n", + " fc = cmap(rescale(\n", + " value, \n", + " np.min(detector_summary[metric]),\n", + " np.max(detector_summary[metric])\n", + " ))\n", + " x, y = zip(*[(corner.x, corner.y) for corner in det.getCorners(cameraGeom.FOCAL_PLANE)])\n", + " ax[0][0].fill(x, y, fc=fc, alpha=alpha)\n", + "\n", + "fig.colorbar(scat, ax=ax[0][0], orientation='vertical', label='Median PSF FWHM (arcsec)')\n", + "ax[0][0].set_xlim(-350, 350)\n", + "ax[0][0].set_ylim(-350, 350)\n", + "ax[0][0].set_xticks([])\n", + "ax[0][0].set_yticks([])\n", + "\n", + "metric = 'psf_e_50'\n", + "\n", + "scat = ax[0][1].scatter(1.e6 + detector_summary['detector'], 1.e6 + detector_summary['detector'], c=detector_summary[metric])\n", + "for det in camera:\n", + " if det.getId() in detector_summary['detector']:\n", + " value = detector_summary.loc[detector_summary['detector'] == det.getId()][metric].values[0]\n", + " fc = cmap(rescale(\n", + " value, \n", + " np.min(detector_summary[metric]),\n", + " np.max(detector_summary[metric])\n", + " ))\n", + " x, y = zip(*[(corner.x, corner.y) for corner in det.getCorners(cameraGeom.FOCAL_PLANE)])\n", + " ax[0][1].fill(x, y, fc=fc, alpha=alpha)\n", + "\n", + "fig.colorbar(scat, ax=ax[0][1], orientation='vertical', label='Median PSF e')\n", + "ax[0][1].set_xlim(-350, 350)\n", + "ax[0][1].set_ylim(-350, 350)\n", + "ax[0][1].set_xticks([])\n", + "ax[0][1].set_yticks([])\n", + "\n", + "metric = 'psf_e1_50'\n", + "\n", + "scat = ax[1][0].scatter(1.e6 + detector_summary['detector'], 1.e6 + detector_summary['detector'], c=detector_summary[metric])\n", + "for det in camera:\n", + " if det.getId() in detector_summary['detector']:\n", + " value = detector_summary.loc[detector_summary['detector'] == det.getId()][metric].values[0]\n", + " fc = cmap(rescale(\n", + " value, \n", + " np.min(detector_summary[metric]),\n", + " np.max(detector_summary[metric])\n", + " ))\n", + " x, y = zip(*[(corner.x, corner.y) for corner in det.getCorners(cameraGeom.FOCAL_PLANE)])\n", + " ax[1][0].fill(x, y, fc=fc, alpha=alpha)\n", + "\n", + "fig.colorbar(scat, ax=ax[1][0], orientation='vertical', label='Median PSF e1')\n", + "ax[1][0].set_xlim(-350, 350)\n", + "ax[1][0].set_ylim(-350, 350)\n", + "ax[1][0].set_xticks([])\n", + "ax[1][0].set_yticks([])\n", + "\n", + "metric = 'psf_e2_50'\n", + "\n", + "scat = ax[1][1].scatter(1.e6 + detector_summary['detector'], 1.e6 + detector_summary['detector'], c=detector_summary[metric])\n", + "for det in camera:\n", + " if det.getId() in detector_summary['detector']:\n", + " value = detector_summary.loc[detector_summary['detector'] == det.getId()][metric].values[0]\n", + " fc = cmap(rescale(\n", + " value, \n", + " np.min(detector_summary[metric]),\n", + " np.max(detector_summary[metric])\n", + " ))\n", + " x, y = zip(*[(corner.x, corner.y) for corner in det.getCorners(cameraGeom.FOCAL_PLANE)])\n", + " ax[1][1].fill(x, y, fc=fc, alpha=alpha)\n", + "\n", + "fig.colorbar(scat, ax=ax[1][1], orientation='vertical', label='Median PSF e2')\n", + "ax[1][1].set_xlim(-350, 350)\n", + "ax[1][1].set_ylim(-350, 350)\n", + "ax[1][1].set_xticks([])\n", + "ax[1][1].set_yticks([])\n", + "\n", + "plt.suptitle('day_obs = %i\\n%s\\n(n = %i)'%(day_obs_select, \n", + " ', '.join(np.unique(ccdvisits['science_program'][selection])),\n", + " np.median(detector_summary['n_visits'])))\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2dba2a73-fedf-4e77-9445-ac89e7a12c66", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(2, 2, figsize=(10, 8), dpi=150)\n", + "\n", + "metric = 'psf_fwhm'\n", + "ax[0][0].errorbar(detector_summary['detector'], \n", + " detector_summary['%s_50'%(metric)],\n", + " yerr=np.array([\n", + " detector_summary['%s_50'%(metric)] - detector_summary['%s_05'%(metric)], \n", + " detector_summary['%s_95'%(metric)] - detector_summary['%s_50'%(metric)],\n", + " ]), \n", + " fmt='_', label='PSF e: 25, 50, 75 percentile')\n", + "ax[0][0].set_xlabel('Detector')\n", + "ax[0][0].set_ylabel('PSF FWHM (arcsec): 25, 50, 75 percentile')\n", + "\n", + "metric = 'psf_e'\n", + "ax[0][1].errorbar(detector_summary['detector'], \n", + " detector_summary['%s_50'%(metric)],\n", + " yerr=np.array([\n", + " detector_summary['%s_50'%(metric)] - detector_summary['%s_05'%(metric)], \n", + " detector_summary['%s_95'%(metric)] - detector_summary['%s_50'%(metric)],\n", + " ]), \n", + " fmt='_', label='PSF e: 25, 50, 75 percentile')\n", + "ax[0][1].set_xlabel('Detector')\n", + "ax[0][1].set_ylabel('PSF e: 25, 50, 75 percentile')\n", + "\n", + "metric = 'psf_e1'\n", + "ax[1][0].errorbar(detector_summary['detector'], \n", + " detector_summary['%s_50'%(metric)],\n", + " yerr=np.array([\n", + " detector_summary['%s_50'%(metric)] - detector_summary['%s_05'%(metric)], \n", + " detector_summary['%s_95'%(metric)] - detector_summary['%s_50'%(metric)],\n", + " ]), \n", + " fmt='_', label='PSF e: 25, 50, 75 percentile')\n", + "ax[1][0].set_xlabel('Detector')\n", + "ax[1][0].set_ylabel('PSF e1: 25, 50, 75 percentile')\n", + "\n", + "metric = 'psf_e2'\n", + "ax[1][1].errorbar(detector_summary['detector'], \n", + " detector_summary['%s_50'%(metric)],\n", + " yerr=np.array([\n", + " detector_summary['%s_50'%(metric)] - detector_summary['%s_05'%(metric)], \n", + " detector_summary['%s_95'%(metric)] - detector_summary['%s_50'%(metric)],\n", + " ]), \n", + " fmt='_', label='PSF e: 25, 50, 75 percentile')\n", + "ax[1][1].set_xlabel('Detector')\n", + "ax[1][1].set_ylabel('PSF e2: 25, 50, 75 percentile')\n", + "\n", + "#plt.suptitle('day_obs = %i\\n(n = %i)'%(day_obs_select, np.median(detector_summary['n_visits'])))\n", + "plt.suptitle('day_obs = %i\\n%s\\n(n = %i)'%(day_obs_select, \n", + " ', '.join(np.unique(ccdvisits['science_program'][selection])),\n", + " np.median(detector_summary['n_visits'])))\n", + "\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "id": "d56045f0-348c-41a2-8e0f-a545ae0b86e8", + "metadata": { + "execution": { + "iopub.execute_input": "2026-04-13T21:46:37.199310Z", + "iopub.status.busy": "2026-04-13T21:46:37.199101Z", + "iopub.status.idle": "2026-04-13T21:46:37.201753Z", + "shell.execute_reply": "2026-04-13T21:46:37.201261Z", + "shell.execute_reply.started": "2026-04-13T21:46:37.199294Z" + } + }, + "source": [ + "# 3rd and 4th order moments" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d610a565-7839-4194-8bdc-747c13fdd93a", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "os.environ[\"no_proxy\"] += \",.consdb\"\n", + "\n", + "from lsst.summit.utils import ConsDbClient\n", + "\n", + "client = ConsDbClient(\"http://consdb-pq.consdb:8080/consdb\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "50c740ea-efd4-4c69-bd63-ec78ff17e177", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import pandas as pd" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2864e45a-acf5-4dd6-b068-85a07bba5a86", + "metadata": {}, + "outputs": [], + "source": [ + "instrument = \"lsstcam\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f8534f84-6cde-41b2-b7ac-52e6bb378f1b", + "metadata": {}, + "outputs": [], + "source": [ + "def fetch(day_obs_min, day_obs_max):\n", + " \n", + " visits_query = f'''\n", + " SELECT \n", + " ccdvisit1_quicklook.psf_area,\n", + " ccdvisit1_quicklook.psf_sigma,\n", + " ccdvisit1_quicklook.psf_ixx,\n", + " ccdvisit1_quicklook.psf_ixy,\n", + " ccdvisit1_quicklook.psf_iyy,\n", + " ccdvisit1_quicklook.coma_1,\n", + " ccdvisit1_quicklook.coma_2,\n", + " ccdvisit1_quicklook.trefoil_1,\n", + " ccdvisit1_quicklook.trefoil_2,\n", + " ccdvisit1_quicklook.e4_1,\n", + " ccdvisit1_quicklook.e4_2,\n", + " ccdvisit1_quicklook.kurtosis, \n", + " ccdvisit1.detector,\n", + " visit1.visit_id,\n", + " visit1.seq_num,\n", + " visit1.band,\n", + " visit1.physical_filter,\n", + " visit1.day_obs,\n", + " visit1.target_name,\n", + " visit1.science_program,\n", + " visit1.observation_reason,\n", + " visit1.airmass,\n", + " visit1.altitude,\n", + " visit1.azimuth,\n", + " visit1.sky_rotation,\n", + " visit1.exp_midpt_mjd,\n", + " visit1.dimm_seeing,\n", + " visit1.s_ra,\n", + " visit1.s_dec,\n", + " visit1_quicklook.psf_sigma_min,\n", + " visit1_quicklook.psf_sigma_median,\n", + " visit1_quicklook.aos_fwhm,\n", + " visit1_quicklook.donut_blur_fwhm,\n", + " visit1_quicklook.ringss_seeing,\n", + " visit1_quicklook.seeing_zenith_500nm_min,\n", + " visit1_quicklook.seeing_zenith_500nm_median,\n", + " visit1_quicklook.physical_rotator_angle\n", + " FROM\n", + " cdb_{instrument}.ccdvisit1_quicklook AS ccdvisit1_quicklook,\n", + " cdb_{instrument}.ccdvisit1 AS ccdvisit1,\n", + " cdb_{instrument}.visit1_quicklook AS visit1_quicklook,\n", + " cdb_{instrument}.visit1 AS visit1 \n", + " WHERE \n", + " ccdvisit1.ccdvisit_id = ccdvisit1_quicklook.ccdvisit_id\n", + " AND ccdvisit1.visit_id = visit1.visit_id \n", + " AND visit1.visit_id = visit1_quicklook.visit_id\n", + " AND ccdvisit1.detector NOT IN (168, 188, 123, 27, 0, 20, 65, 161) -- vignetted\n", + " AND ccdvisit1.detector NOT IN (191, 192, 195, 196, 199, 200, 203, 204) -- corner wavefront sensors\n", + " AND visit1.airmass > 0\n", + " AND visit1.day_obs >= {day_obs_min} AND visit1.day_obs <= {day_obs_max}\n", + " -- AND visit1.science_program in ('BLOCK-365', 'BLOCK-407');\n", + " AND visit1.science_program in ('BLOCK-365', 'BLOCK-407', 'BLOCK-408', 'BLOCK-416', 'BLOCK-417', 'BLOCK-419', 'BLOCK-421') \n", + " -- , 'BLOCK-T698')\n", + " -- AND visit1.science_program in ('BLOCK-T539', 'BLOCK-T698', 'BLOCK-T641') -- initial alignment, closed loop stability\n", + " AND visit1.physical_filter in ('u_24', 'g_6', 'r_57', 'i_39', 'z_20', 'y_10');\n", + " '''\n", + " \n", + " ccdvisits = client.query(visits_query).to_pandas()\n", + "\n", + " pixel_scale = 0.2\n", + " sig2fwhm = 2 * np.sqrt(2 * np.log(2))\n", + " ccdvisits[\"psf_fwhm\"] = ccdvisits[\"psf_sigma\"] * sig2fwhm * pixel_scale\n", + " ccdvisits[\"psf_fwhm\"] = pd.to_numeric(ccdvisits[\"psf_fwhm\"], errors=\"coerce\")\n", + " ccdvisits[\"donut_blur_fwhm\"] = pd.to_numeric(ccdvisits[\"donut_blur_fwhm\"], errors=\"coerce\")\n", + " ccdvisits[\"aos_fwhm\"] = pd.to_numeric(ccdvisits[\"aos_fwhm\"], errors=\"coerce\")\n", + "\n", + " ccdvisits[\"psf_area\"] = pd.to_numeric(ccdvisits[\"psf_area\"], errors=\"coerce\")\n", + " ccdvisits[\"psf_fwhm_area\"] = 0.663 * pixel_scale * np.sqrt(ccdvisits[\"psf_area\"])\n", + "\n", + " #airmass_seeing_correction = np.array([getAirmassSeeingCorrection(airmass) for airmass in ccdvisits[\"airmass\"]])\n", + " #bandpass_seeing_correction = np.array([getBandpassSeeingCorrection(band) for band in ccdvisits[\"physical_filter\"]])\n", + " #ccdvisits[\"psf_fwhm_zenith_500nm\"] = ccdvisits[\"psf_fwhm\"] * airmass_seeing_correction * bandpass_seeing_correction\n", + "\n", + " #ccdvisits[\"psf_e1\"] = (ccdvisits[\"psf_ixx\"]**2 - ccdvisits[\"psf_iyy\"]**2) / (ccdvisits[\"psf_ixx\"]**2 + ccdvisits[\"psf_iyy\"]**2)\n", + " #ccdvisits[\"psf_e2\"] = (2. * ccdvisits[\"psf_ixy\"]**2) / (ccdvisits[\"psf_ixx\"]**2 + ccdvisits[\"psf_iyy\"]**2)\n", + " ccdvisits[\"psf_e1\"] = (ccdvisits[\"psf_ixx\"] - ccdvisits[\"psf_iyy\"]) / (ccdvisits[\"psf_ixx\"] + ccdvisits[\"psf_iyy\"])\n", + " ccdvisits[\"psf_e2\"] = (2. * ccdvisits[\"psf_ixy\"]) / (ccdvisits[\"psf_ixx\"] + ccdvisits[\"psf_iyy\"])\n", + " ccdvisits[\"psf_e\"] = np.sqrt(np.array(ccdvisits[\"psf_e1\"]**2 + ccdvisits[\"psf_e2\"]**2, dtype=float))\n", + " ccdvisits[\"psf_e1_unnormalized\"] = (ccdvisits[\"psf_ixx\"] - ccdvisits[\"psf_iyy\"])\n", + " ccdvisits[\"psf_e2_unnormalized\"] = (2. * ccdvisits[\"psf_ixy\"])\n", + " ccdvisits[\"psf_e_unnormalized\"] = np.sqrt(np.array(ccdvisits[\"psf_e1_unnormalized\"]**2 + ccdvisits[\"psf_e2_unnormalized\"]**2, dtype=float))\n", + "\n", + " ccdvisits[\"psf_e\"] = pd.to_numeric(ccdvisits[\"psf_e\"], errors=\"coerce\")\n", + " ccdvisits[\"psf_e1\"] = pd.to_numeric(ccdvisits[\"psf_e1\"], errors=\"coerce\")\n", + " ccdvisits[\"psf_e2\"] = pd.to_numeric(ccdvisits[\"psf_e2\"], errors=\"coerce\")\n", + "\n", + " ccdvisits[\"coma_1\"] = pd.to_numeric(ccdvisits[\"coma_1\"], errors=\"coerce\")\n", + " ccdvisits[\"coma_2\"] = pd.to_numeric(ccdvisits[\"coma_2\"], errors=\"coerce\")\n", + " ccdvisits[\"trefoil_1\"] = pd.to_numeric(ccdvisits[\"trefoil_1\"], errors=\"coerce\")\n", + " ccdvisits[\"trefoil_2\"] = pd.to_numeric(ccdvisits[\"trefoil_2\"], errors=\"coerce\")\n", + "\n", + " return ccdvisits" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d50ec475-7294-41ca-a84a-2b627d7c3786", + "metadata": {}, + "outputs": [], + "source": [ + "day_obs_min = 20260401\n", + "day_obs_max = 20260421\n", + "ccdvisits = fetch(day_obs_min, day_obs_max)\n", + "print(len(ccdvisits))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ec70f9c2-db67-45b8-a151-a9d65fbde011", + "metadata": {}, + "outputs": [], + "source": [ + "ccdvisits.columns" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1e72a161-5658-4ac4-afa7-46bfa3e3977c", + "metadata": {}, + "outputs": [], + "source": [ + "np.any(ccdvisits['trefoil_1'].notna())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "23485b66-fb2f-4c60-af82-58851f49766f", + "metadata": {}, + "outputs": [], + "source": [ + "ccdvisits[['visit_id', 'coma_1', 'coma_2', 'trefoil_1', 'trefoil_2', 'e4_1', 'e4_2', 'kurtosis']]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "34570f28-5639-4cf5-8d5b-fab0e5d1df22", + "metadata": {}, + "outputs": [], + "source": [ + "# https://rubin-obs.slack.com/archives/C0AUAM3UPB2/p1776637827811059?thread_ts=1776631306.523939&cid=C0AUAM3UPB2\n", + "# moments_score = 3*(coma_1**2 + coma_2**2) + (trefoil_1**2 + trefoil_2**2)\n", + "# Estimate of non-Gaussian optical power from Gaussian weighted 3rd order PSF star moments\n", + "\n", + "ccdvisits['coma'] = np.sqrt(ccdvisits['coma_1']**2 + ccdvisits['coma_2']**2)\n", + "ccdvisits['trefoil'] = np.sqrt(ccdvisits['trefoil_1']**2 + ccdvisits['trefoil_2']**2)\n", + "ccdvisits['moment_score'] = 3 * ccdvisits['coma']**2 + ccdvisits['trefoil']**2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "794b500f-9d07-498e-9fc6-e488ff62438c", + "metadata": {}, + "outputs": [], + "source": [ + "print(0.005 / 0.001)\n", + "print(0.15 / 0.035)\n", + "print(0.2 / 0.055)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8c853de5-cde9-4bb4-a619-d64de8dc49fb", + "metadata": {}, + "outputs": [], + "source": [ + "log = True\n", + "bins = np.linspace(0, 0.2, 1001)\n", + "\n", + "plt.figure()\n", + "plt.hist(ccdvisits['moment_score'], bins=bins, cumulative=-1, density=True, histtype='step')\n", + "\n", + "thresholds = 4. * np.array([0.002, 0.005, 0.02])\n", + "print(thresholds)\n", + "for threshold in thresholds:\n", + " plt.axvline(threshold, color='0.5', ls='--', lw=1)\n", + "\n", + " fraction = np.sum(ccdvisits['moment_score'] > threshold) / float(len(ccdvisits['moment_score']))\n", + " \n", + " plt.axhline(fraction, color='0.5', ls='--', lw=1)\n", + "\n", + "plt.xlim(bins[0], 0.2)\n", + "\n", + "if log:\n", + " plt.ylim(0.0001, 1.)\n", + " plt.yscale('log')\n", + "else:\n", + " plt.ylim(0., 1.)\n", + "plt.xlabel('Moment Score')\n", + "plt.ylabel('1 - CDF')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7b77f4a1-3231-4468-ad2c-105411ea6219", + "metadata": {}, + "outputs": [], + "source": [ + "selection = (ccdvisits['donut_blur_fwhm'] < 0.8)\n", + "\n", + "groups = ccdvisits[selection].groupby('day_obs')\n", + "day_obs_summary = pd.DataFrame({\n", + " 'day_obs': groups['day_obs'].first(),\n", + " 'n_ccdvisits': groups['day_obs'].count(),\n", + " 'moment_score_low': groups['moment_score'].agg(lambda x: (x <= 4. * 0.002).mean()),\n", + " 'moment_score_002': groups['moment_score'].agg(lambda x: (x > 4. * 0.002).mean()),\n", + " 'moment_score_005': groups['moment_score'].agg(lambda x: (x > 4. * 0.005).mean()),\n", + " 'moment_score_02': groups['moment_score'].agg(lambda x: (x > 4. * 0.02).mean()),\n", + "})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8e4b6e36-fb49-43ea-835f-3882a5d5b1ca", + "metadata": {}, + "outputs": [], + "source": [ + "day_obs_summary" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9040dd81-3875-439f-8e68-5ea8f445d119", + "metadata": {}, + "outputs": [], + "source": [ + "from astropy.time import Time\n", + "\n", + "def dayObsToTime(day_obs):\n", + " year = str(day_obs)[0:4]\n", + " month = str(day_obs)[4:6]\n", + " day = str(day_obs)[6:8]\n", + " time = Time('%s-%s-%s 12:00:00'%(year, month, day), format='iso')\n", + " return time\n", + "\n", + "#selection = np.tile(True, len(day_obs_summary))\n", + "selection = (day_obs_summary['n_ccdvisits'] >= (21 * 10))\n", + "\n", + "xticks = np.array([dayObsToTime(_).mjd for _ in day_obs_summary['day_obs'][selection]])\n", + "\n", + "print(xticks)\n", + "\n", + "fig, ax = plt.subplots(figsize=(15, 6), dpi=150)\n", + "a = np.linspace(0., 1., 10)\n", + "\n", + "ax.scatter(xticks + 0., day_obs_summary['moment_score_low'][selection], c='tab:green', label='Moment Score < 0.008')\n", + "#ax.scatter(xticks + 0., 1. - day_obs_summary['shapelets_score_002'][selection])\n", + "ax.scatter(xticks + 0., 1. - day_obs_summary['moment_score_005'][selection], c='tab:orange', label='Moment Score < 0.02')\n", + "#ax.scatter(xticks + 0., day_obs_summary['shapelets_score_02'][selection])\n", + "\n", + "#first_legend = ax.legend(handles=[variation, aos, sys], loc='upper left', title='Estimated System Contribution')\n", + "#second_legend = ax.legend(handles=[psf, atm, ellipticity], loc='upper right', title='Measured Quantities')\n", + "#ax.add_artist(first_legend)\n", + "ax.legend(loc='lower right')\n", + "\n", + "ax.set_xticks(xticks, day_obs_summary['day_obs'][selection], rotation=90, fontsize=7.5)\n", + "ax.set_ylim(0., 1.0)\n", + "\n", + "#ax.axhline(0.45, c='red', ls=':', zorder=0)\n", + "#ax.axhline(0.04, c='orange', ls=':', zorder=0)\n", + "\n", + "#ax.set_ylabel('Image Quality Metric (arcsec)')\n", + "ax.set_ylabel('Fraction')\n", + "#plt.title('10th, 50th, 90th Percentiles of Per-visit Quantities for Each Night')\n", + "plt.grid()\n", + "plt.title('Donut Blur < 0.8')\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e0abf2f5-3f88-4cf0-b592-7b5d2124d9e8", + "metadata": {}, + "outputs": [], + "source": [ + "log = False\n", + "bins = np.linspace(0, 0.2, 1001)\n", + "\n", + "selection = (ccdvisits['donut_blur_fwhm'] < 0.8)\n", + "\n", + "plt.figure()\n", + "#plt.hist(merged_df['shapelets_score'], bins=bins, cumulative=-1, density=True, color='black', histtype='step')\n", + "plt.hist(ccdvisits['moment_score'][selection], bins=bins, cumulative=-1, density=True, histtype='step', label='Donut Blur < 0.8')\n", + "plt.hist(ccdvisits['moment_score'][~selection], bins=bins, cumulative=-1, density=True, histtype='step', label='Donut Blur > 0.8')\n", + "\n", + "thresholds = 4. * np.array([0.002, 0.005, 0.02])\n", + "print(thresholds)\n", + "for threshold in thresholds:\n", + " plt.axvline(threshold, color='0.5', ls='--', lw=1)\n", + "\n", + " fraction = np.sum(ccdvisits['moment_score'] > threshold) / float(len(ccdvisits['moment_score']))\n", + " \n", + " #plt.axhline(fraction, color='0.5', ls='--', lw=1)\n", + "\n", + "plt.xlim(bins[0], 0.2)\n", + "\n", + "plt.axvspan(0., 4. * 0.002, color='tab:blue', alpha=0.1)\n", + "plt.axvspan(4. * 0.002, 4. * 0.005, color='tab:green', alpha=0.1)\n", + "plt.axvspan(4. * 0.005, 4. * 0.02, color='yellow', alpha=0.1)\n", + "plt.axvspan(4. * 0.02, 4. * 0.05, color='tab:red', alpha=0.1)\n", + "\n", + "if log:\n", + " plt.ylim(0.0001, 1.)\n", + " plt.yscale('log')\n", + "else:\n", + " plt.ylim(0., 1.)\n", + "plt.legend()\n", + "plt.xlabel('Moment Score')\n", + "plt.ylabel('1 - CDF')\n", + "plt.title('%s - %s'%(day_obs_min, day_obs_max))\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "6c952d22-e412-4441-9348-6b5b436ceae5", + "metadata": { + "execution": { + "iopub.execute_input": "2026-04-17T18:44:34.867122Z", + "iopub.status.busy": "2026-04-17T18:44:34.866892Z", + "iopub.status.idle": "2026-04-17T18:44:34.869564Z", + "shell.execute_reply": "2026-04-17T18:44:34.868960Z", + "shell.execute_reply.started": "2026-04-17T18:44:34.867107Z" + }, + "jp-MarkdownHeadingCollapsed": true + }, + "source": [ + "# Look at Shapelets" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f06fae00-c40f-4fae-a10d-58a36e0768a5", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import pandas as pd\n", + "\n", + "import matplotlib.pyplot as plt\n", + "\n", + "%matplotlib widget" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b564a5bc-ed9f-4f30-99a4-9410b1ebcc94", + "metadata": {}, + "outputs": [], + "source": [ + "infile = '/sdf/data/rubin/user/ztq1996/psf-rubin/psf_zernike/data/shapelets_score_all_visits.pq'\n", + "\n", + "df_shapelet = pd.read_parquet(infile)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9bb6b2a6-a70c-4826-b195-7db5f5527ef7", + "metadata": {}, + "outputs": [], + "source": [ + "df_shapelet" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fe8c1f99-f8fd-43c5-a163-b793c44994b0", + "metadata": {}, + "outputs": [], + "source": [ + "df_shapelet.columns" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12ae6aee-2965-4a03-a7f0-c1ec7efe9ae2", + "metadata": {}, + "outputs": [], + "source": [ + "df_shapelet['visit_id']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "034958f4-3752-4c7c-b30a-98075a0b208b", + "metadata": {}, + "outputs": [], + "source": [ + "log = True\n", + "bins = np.linspace(0, 0.2, 1001)\n", + "\n", + "plt.figure()\n", + "plt.hist(df_shaplet['shapelets_score'], bins=bins, cumulative=-1, density=True, histtype='step')\n", + "\n", + "for threshold in [0.002, 0.005, 0.02]:\n", + " plt.axvline(threshold, color='0.5', ls='--', lw=1)\n", + "\n", + " fraction = np.sum(df_shaplet['shapelets_score'] > threshold) / float(len(df_shaplet['shapelets_score']))\n", + " \n", + " plt.axhline(fraction, color='0.5', ls='--', lw=1)\n", + "\n", + "plt.xlim(bins[0], 0.04)\n", + "\n", + "if log:\n", + " plt.ylim(0.0001, 1.)\n", + " plt.yscale('log')\n", + "else:\n", + " plt.ylim(0., 1.)\n", + "plt.xlabel('Shapelet Score')\n", + "plt.ylabel('1 - CDF')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bf461fbf-858e-4e02-8a8c-74466e4eb846", + "metadata": {}, + "outputs": [], + "source": [ + "def fetch(day_obs_min, day_obs_max):\n", + " \n", + " visits_query = f'''\n", + " SELECT \n", + " ccdvisit1.detector,\n", + " visit1.visit_id,\n", + " visit1.seq_num,\n", + " visit1.band,\n", + " visit1.physical_filter,\n", + " visit1.day_obs,\n", + " visit1.target_name,\n", + " visit1.science_program,\n", + " visit1.observation_reason,\n", + " visit1.airmass,\n", + " visit1.altitude,\n", + " visit1.azimuth,\n", + " visit1.sky_rotation,\n", + " visit1.exp_midpt_mjd,\n", + " visit1.dimm_seeing,\n", + " visit1.s_ra,\n", + " visit1.s_dec,\n", + " visit1_quicklook.donut_blur_fwhm\n", + " FROM\n", + " cdb_{instrument}.ccdvisit1 AS ccdvisit1,\n", + " cdb_{instrument}.visit1 AS visit1, \n", + " cdb_{instrument}.visit1_quicklook AS visit1_quicklook\n", + " WHERE \n", + " ccdvisit1.visit_id = visit1.visit_id\n", + " AND visit1.visit_id = visit1_quicklook.visit_id\n", + " AND ccdvisit1.detector NOT IN (168, 188, 123, 27, 0, 20, 65, 161) -- vignetted\n", + " AND ccdvisit1.detector NOT IN (191, 192, 195, 196, 199, 200, 203, 204) -- corner wavefront sensors\n", + " AND visit1.airmass > 0\n", + " AND visit1.day_obs >= {day_obs_min} AND visit1.day_obs <= {day_obs_max}\n", + " -- AND visit1.science_program in ('BLOCK-365', 'BLOCK-407');\n", + " AND visit1.science_program in ('BLOCK-365', 'BLOCK-407', 'BLOCK-408', 'BLOCK-416', 'BLOCK-417', 'BLOCK-419', 'BLOCK-421') \n", + " -- , 'BLOCK-T698')\n", + " -- AND visit1.science_program in ('BLOCK-T539', 'BLOCK-T698', 'BLOCK-T641') -- initial alignment, closed loop stability\n", + " AND visit1.physical_filter in ('u_24', 'g_6', 'r_57', 'i_39', 'z_20', 'y_10');\n", + " '''\n", + " \n", + " ccdvisits = client.query(visits_query).to_pandas()\n", + "\n", + " return ccdvisits" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7682f85d-a21a-4abb-bf30-d0a51a11a25b", + "metadata": {}, + "outputs": [], + "source": [ + "#day_obs_min = 20250620 # Start of SV surveys\n", + "#day_obs_max = 20250921 # Last night of SV surveys\n", + "day_obs_min = 20251026\n", + "day_obs_max = 20260415\n", + "ccdvisits = fetch(day_obs_min, day_obs_max)\n", + "print(len(ccdvisits))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6652e38d-1624-414c-8db0-d14362385ad4", + "metadata": {}, + "outputs": [], + "source": [ + "merged_df = pd.merge(ccdvisits, df_shapelet, \n", + " left_on=['visit_id', 'detector'], \n", + " right_on=['visit_id', 'detector_id'])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1fe39c73-0b3e-465d-933c-71c74438638e", + "metadata": {}, + "outputs": [], + "source": [ + "merged_df" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2b71fb2a-39d2-400c-93e3-6e0c955bd99f", + "metadata": { + "jupyter": { + "source_hidden": true + } + }, + "outputs": [], + "source": [ + "#groups = ccdvisits.groupby('visit_id')\n", + "groups = merged_df.groupby('visit_id')\n", + "visits_summary = pd.DataFrame({\n", + " 'day_obs': groups['day_obs'].first(),\n", + " 'target_name': groups['target_name'].first(),\n", + " 'science_program': groups['science_program'].first(),\n", + " 'observation_reason': groups['observation_reason'].first(),\n", + " 'seq_num': groups['seq_num'].median(),\n", + " 'exp_midpt_mjd': groups['exp_midpt_mjd'].median(),\n", + " 'donut_blur_fwhm': groups['donut_blur_fwhm'].median(),\n", + " 'aos_fwhm': groups['aos_fwhm'].median(),\n", + " 'donut_blur_atm_fwhm': groups['donut_blur_atm_fwhm'].median(),\n", + " 'aos_cam_fwhm': groups['aos_cam_fwhm'].median(),\n", + " 'physical_rotator_angle': groups['physical_rotator_angle'].median(),\n", + " 'altitude': groups['altitude'].median(),\n", + " 'airmass': groups['airmass'].median(),\n", + " 'psf_fwhm_05': groups['psf_fwhm'].quantile(0.05),\n", + " 'psf_fwhm_50': groups['psf_fwhm'].quantile(0.50),\n", + " 'psf_fwhm_95': groups['psf_fwhm'].quantile(0.95),\n", + " 'psf_fwhm_zenith_500nm_50': groups['psf_fwhm_zenith_500nm'].quantile(0.50),\n", + " 'psf_fwhm_area_50': groups['psf_fwhm_area'].quantile(0.50),\n", + " 'psf_e_05': groups['psf_e'].quantile(0.05),\n", + " 'psf_e_50': groups['psf_e'].quantile(0.50),\n", + " 'psf_e_95': groups['psf_e'].quantile(0.95),\n", + " 'psf_e1_05': groups['psf_e1'].quantile(0.05),\n", + " 'psf_e1_50': groups['psf_e1'].quantile(0.50),\n", + " 'psf_e1_95': groups['psf_e1'].quantile(0.95),\n", + " 'psf_e2_05': groups['psf_e2'].quantile(0.05),\n", + " 'psf_e2_50': groups['psf_e2'].quantile(0.50),\n", + " 'psf_e2_95': groups['psf_e2'].quantile(0.95),\n", + " 'psf_e_unnormalized_05': groups['psf_e_unnormalized'].quantile(0.05),\n", + " 'psf_e_unnormalized_50': groups['psf_e_unnormalized'].quantile(0.50),\n", + " 'psf_e_unnormalized_95': groups['psf_e_unnormalized'].quantile(0.95),\n", + " 'band': groups['band'].first(),\n", + " 'shapelets_score_05': groups['shapelets_score'].quantile(0.05),\n", + " 'shapelets_score_50': groups['shapelets_score'].quantile(0.50),\n", + " 'shapelets_score_95': groups['shapelets_score'].quantile(0.95),\n", + " 'n_ccdvisits': groups['detector'].count(),\n", + "})\n", + "visits_summary['psf_fwhm_95_05'] = np.sqrt(visits_summary['psf_fwhm_95']**2 - visits_summary['psf_fwhm_05']**2)\n", + "visits_summary['sys_50'] = np.sqrt(visits_summary['psf_fwhm_50']**2 + CAM_FWHM**2 - visits_summary['donut_blur_fwhm']**2)\n", + "visits_summary['psf_fwhm_model'] = np.sqrt(visits_summary['aos_fwhm']**2 + visits_summary['donut_blur_fwhm']**2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "05e98015-931d-40d5-8ded-d7e5d508beaf", + "metadata": {}, + "outputs": [], + "source": [ + "groups = merged_df.groupby('visit_id')\n", + "visits_summary = pd.DataFrame({\n", + " 'n_ccdvisits': groups['detector'].count(),\n", + "})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7a4fbfe6-4b24-4614-8c0c-dcef2e905d63", + "metadata": {}, + "outputs": [], + "source": [ + "#np.max(visits_summary)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dc053d03-3b8f-4cc0-9e72-500bce7dc3f1", + "metadata": {}, + "outputs": [], + "source": [ + "#print(merged_df['detector'][merged_df['visit_id'] == 2025112400345])\n", + "#print(merged_df['detector'][merged_df['visit_id'] == 2026040900194])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cfdeb31c-0650-4269-a5fc-1893597acb98", + "metadata": { + "jupyter": { + "source_hidden": true + } + }, + "outputs": [], + "source": [ + "day_obs = 20260414\n", + "selection = (visits_summary['day_obs'] == day_obs)\n", + "\n", + "print(np.sum(selection))\n", + "\n", + "fig, ax = plt.subplots(3, sharex=True, figsize=(16, 10))\n", + "\n", + "ax[0].errorbar(visits_summary['seq_num'][selection], \n", + " visits_summary['psf_fwhm_50'][selection],\n", + " yerr=np.array([\n", + " visits_summary['psf_fwhm_50'][selection] - visits_summary['psf_fwhm_05'][selection], \n", + " visits_summary['psf_fwhm_95'][selection] - visits_summary['psf_fwhm_50'][selection],\n", + " ]), \n", + " fmt='_', label='PSF FWHM: 5, 50, 95 percentile')\n", + "ax[0].errorbar(visits_summary['seq_num'][selection], \n", + " visits_summary['psf_e_50'][selection],\n", + " yerr=np.array([\n", + " visits_summary['psf_e_50'][selection] - visits_summary['psf_e_05'][selection], \n", + " visits_summary['psf_e_95'][selection] - visits_summary['psf_e_50'][selection],\n", + " ]), \n", + " fmt='_', label='PSF e: 5, 50, 95 percentile')\n", + "\n", + "shapelests_scale = 100.\n", + "ax[0].errorbar(visits_summary['seq_num'][selection], \n", + " shapelests_scale * visits_summary['shapelets_score_50'][selection],\n", + " yerr=np.array([\n", + " shapelests_scale * (visits_summary['shapelets_score_50'][selection] - visits_summary['shapelets_score_05'][selection]), \n", + " shapelests_scale * (visits_summary['shapelets_score_95'][selection] - visits_summary['shapelets_score_50'][selection]),\n", + " ]), \n", + " fmt='_', label='Shapelets Score: 5, 50, 95 percentile', c='magenta')\n", + "\n", + "\n", + "#ax[0].errorbar(visits_summary['seq_num'][selection] + 0.2, \n", + "# 0.2 * visits_summary['psf_e_unnormalized_50'][selection],\n", + "# yerr=np.array([\n", + "# 0.2 * (visits_summary['psf_e_unnormalized_50'][selection] - visits_summary['psf_e_unnormalized_05'][selection]), \n", + "# 0.2 * (visits_summary['psf_e_unnormalized_95'][selection] - visits_summary['psf_e_unnormalized_50'][selection]),\n", + "# ]), \n", + "# fmt='_', label='PSF e (unnormalized): 5, 50, 95 percentile', c='magenta')\n", + "#ax[0].scatter(visits_summary['seq_num'], visits_summary['psf_fwhm_95'])\n", + "#ax[0].scatter(visits_summary['seq_num'], visits_summary['psf_fwhm_50'])\n", + "#ax[0].scatter(visits_summary['seq_num'], visits_summary['psf_fwhm_05'])\n", + "ax[0].scatter(visits_summary['seq_num'][selection], visits_summary['psf_fwhm_95_05'][selection], c='red', marker='x', label='sqrt(fwhm_95^2 - fwhm_5^2)')\n", + "ax[0].scatter(visits_summary['seq_num'][selection], visits_summary['donut_blur_fwhm'][selection], c='black', marker='o', label='Donut Blur')\n", + "ax[0].scatter(visits_summary['seq_num'][selection], visits_summary['psf_fwhm_model'][selection], c='0.8', marker='o', label='Model')\n", + "ax[0].scatter(visits_summary['seq_num'][selection], visits_summary['aos_fwhm'][selection], c='green', marker='+', label='AOS Residual')\n", + "#ax[0].scatter(visits_summary['seq_num'][selection], np.fabs(visits_summary[\"z4_fwhm\"][selection] + 0.25), c='magenta', marker='o', label='z4_fwhm')\n", + "ax[0].set_ylabel('Delivered Image Quality Metric\\n(arcsec)')\n", + "#ax[0].scatter(visit_seq_num, np.sqrt(visit_donut_blur_fwhm**2 + visit_psf_fwhm_95_05**2), c='green', marker='+', label='AOS Residual')\n", + "ax[0].legend(loc='upper right')\n", + "#ax[0].set_ylim(0., 2.0)\n", + "ax[0].set_ylim(0., 1.4)\n", + "ax[0].grid(True)\n", + "\n", + "for band in bands:\n", + " band_selection = selection & (visits_summary[\"band\"] == band)\n", + " ax[1].scatter(visits_summary['seq_num'][band_selection], visits_summary['altitude'][band_selection], marker='.', s=50, c=colors[band], label=band)\n", + "#ax[1].scatter(visits_summary['seq_num'], visits_summary['altitude'], marker='.', s=10, c='black')\n", + "ax[1].set_ylabel('Elevation\\n(deg)')\n", + "ax[1].grid(True)\n", + "ax[1].legend()\n", + "\n", + "for band in bands:\n", + " band_selection = selection & (visits_summary[\"band\"] == band)\n", + " ax[2].scatter(visits_summary['seq_num'][band_selection], visits_summary['physical_rotator_angle'][band_selection], marker='.', s=50, c=colors[band], label=band)\n", + "#ax[2].scatter(visits_summary['seq_num'], visits_summary['physical_rotator_angle'], marker='.', s=10, c='black')\n", + "ax[2].set_ylabel('Phys Rot Angle\\n(deg)')\n", + "ax[2].grid(True)\n", + "\n", + "ax[-1].set_xlabel('seq_num')\n", + "\n", + "plt.suptitle('day_obs = %s'%(day_obs))\n", + "\n", + "#plt.suptitle('day_obs = %s'%(day_obs_min))\n", + "plt.subplots_adjust(hspace=0)\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a01db4d1-578b-4c7e-886e-bd9df1154e84", + "metadata": {}, + "outputs": [], + "source": [ + "np.sum(merged_df['visit_id'] == 2026040900194)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ed65fe76-6f27-4f3a-a5fa-56a086fb1dd6", + "metadata": {}, + "outputs": [], + "source": [ + "# [0.002, 0.005, 0.02]\n", + "\n", + "selection = (merged_df['donut_blur_fwhm'] < 0.8)\n", + "\n", + "groups = merged_df[selection].groupby('day_obs')\n", + "day_obs_summary = pd.DataFrame({\n", + " 'day_obs': groups['day_obs'].first(),\n", + " 'n_ccdvisits': groups['day_obs'].count(),\n", + " 'shapelets_score_low': groups['shapelets_score'].agg(lambda x: (x <= 0.002).mean()),\n", + " 'shapelets_score_002': groups['shapelets_score'].agg(lambda x: (x > 0.002).mean()),\n", + " 'shapelets_score_005': groups['shapelets_score'].agg(lambda x: (x > 0.005).mean()),\n", + " 'shapelets_score_02': groups['shapelets_score'].agg(lambda x: (x > 0.02).mean()),\n", + "})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f62d7958-5af6-4288-9f35-adf21a3bd57e", + "metadata": {}, + "outputs": [], + "source": [ + "day_obs_summary" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d44f8435-8d50-4f7d-86bc-f83250462061", + "metadata": {}, + "outputs": [], + "source": [ + "day_obs_summary['n_ccdvisits'].values" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "64b7e58d-97ce-4d28-85b8-449a17e16331", + "metadata": {}, + "outputs": [], + "source": [ + "from astropy.time import Time\n", + "\n", + "def dayObsToTime(day_obs):\n", + " year = str(day_obs)[0:4]\n", + " month = str(day_obs)[4:6]\n", + " day = str(day_obs)[6:8]\n", + " time = Time('%s-%s-%s 12:00:00'%(year, month, day), format='iso')\n", + " return time\n", + "\n", + "#selection = np.tile(True, len(day_obs_summary))\n", + "selection = (day_obs_summary['n_ccdvisits'] >= (21 * 10))\n", + "\n", + "xticks = np.array([dayObsToTime(_).mjd for _ in day_obs_summary['day_obs'][selection]])\n", + "\n", + "print(xticks)\n", + "\n", + "fig, ax = plt.subplots(figsize=(15, 6), dpi=150)\n", + "a = np.linspace(0., 1., 10)\n", + "\n", + "ax.scatter(xticks + 0., day_obs_summary['shapelets_score_low'][selection], c='tab:green', label='Shapelets Score < 0.002')\n", + "#ax.scatter(xticks + 0., 1. - day_obs_summary['shapelets_score_002'][selection])\n", + "ax.scatter(xticks + 0., 1. - day_obs_summary['shapelets_score_005'][selection], c='tab:orange', label='Shapelets Score < 0.005')\n", + "#ax.scatter(xticks + 0., day_obs_summary['shapelets_score_02'][selection])\n", + "\n", + "#first_legend = ax.legend(handles=[variation, aos, sys], loc='upper left', title='Estimated System Contribution')\n", + "#second_legend = ax.legend(handles=[psf, atm, ellipticity], loc='upper right', title='Measured Quantities')\n", + "#ax.add_artist(first_legend)\n", + "ax.legend(loc='lower right')\n", + "\n", + "ax.set_xticks(xticks, day_obs_summary['day_obs'][selection], rotation=90, fontsize=7.5)\n", + "ax.set_ylim(0., 1.0)\n", + "\n", + "#ax.axhline(0.45, c='red', ls=':', zorder=0)\n", + "#ax.axhline(0.04, c='orange', ls=':', zorder=0)\n", + "\n", + "#ax.set_ylabel('Image Quality Metric (arcsec)')\n", + "ax.set_ylabel('Fraction')\n", + "#plt.title('10th, 50th, 90th Percentiles of Per-visit Quantities for Each Night')\n", + "plt.grid()\n", + "plt.title('Donut Blur < 0.8')\n", + "fig.tight_layout()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e79a0867-0ec1-4142-ab5f-10b65a147152", + "metadata": {}, + "outputs": [], + "source": [ + "log = False\n", + "bins = np.linspace(0, 0.2, 1001)\n", + "\n", + "selection = (merged_df['donut_blur_fwhm'] < 0.8)\n", + "\n", + "plt.figure()\n", + "#plt.hist(merged_df['shapelets_score'], bins=bins, cumulative=-1, density=True, color='black', histtype='step')\n", + "plt.hist(merged_df['shapelets_score'][selection], bins=bins, cumulative=-1, density=True, histtype='step', label='Donut Blur < 0.8')\n", + "plt.hist(merged_df['shapelets_score'][~selection], bins=bins, cumulative=-1, density=True, histtype='step', label='Donut Blur > 0.8')\n", + "\n", + "for threshold in [0.002, 0.005, 0.02]:\n", + " plt.axvline(threshold, color='0.5', ls='--', lw=1)\n", + "\n", + " fraction = np.sum(merged_df['shapelets_score'] > threshold) / float(len(merged_df['shapelets_score']))\n", + " \n", + " #plt.axhline(fraction, color='0.5', ls='--', lw=1)\n", + "\n", + "plt.xlim(bins[0], 0.04)\n", + "\n", + "plt.axvspan(0., 0.002, color='tab:blue', alpha=0.1)\n", + "plt.axvspan(0.002, 0.005, color='tab:green', alpha=0.1)\n", + "plt.axvspan(0.005, 0.02, color='yellow', alpha=0.1)\n", + "plt.axvspan(0.02, 0.04, color='tab:red', alpha=0.1)\n", + "\n", + "if log:\n", + " plt.ylim(0.0001, 1.)\n", + " plt.yscale('log')\n", + "else:\n", + " plt.ylim(0., 1.)\n", + "plt.legend()\n", + "plt.xlabel('Shapelet Score')\n", + "plt.ylabel('1 - CDF')\n", + "plt.title('%s - %s'%(day_obs_min, day_obs_max))\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "40b26035-5e5a-4d03-8735-8dbeef497b1b", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "LSST", + "language": "python", + "name": "lsst" + }, + "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.13.9" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/notebooks/survey_performance.ipynb b/notebooks/survey_performance.ipynb deleted file mode 100644 index 0b603eb..0000000 --- a/notebooks/survey_performance.ipynb +++ /dev/null @@ -1,41 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Current Estimate of Survey Performance for {{params.instrument}} on {{ params.date }}" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Parameters. Set defaults here.\n", - "\n", - "date = \"2024-12-11\"\n", - "instrument = \"LSSTComCam\"" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import numpy as np\n", - "import matplotlib.pyplot as plt\n", - "from lsst_efd_client import EfdClient" - ] - } - ], - "metadata": { - "language_info": { - "name": "python" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/notebooks/survey_performance.yaml b/notebooks/survey_performance.yaml deleted file mode 100644 index ab98098..0000000 --- a/notebooks/survey_performance.yaml +++ /dev/null @@ -1,18 +0,0 @@ -title: Survey Performance Estimate -description: Current estimate of the Survey Performance -authors: - - name: Leanne Guy, Željko Ivezić - slack: leannep -parameters: - date: - type: string - description: Date of summary (YYYY-MM-DD) - default: "2024-09-04" - instrument: - type: string - description: Instrument ("LSSTCam", "LSSTComCam") - default: "LSSTCam" - telescope: - type: string - description: The telescope to compute the performance for - default: "Simonyi"