Fix uncertainty error

desy-multimessenger · Jun 30, 2022 · 4d9a3f2 · 4d9a3f2
1 parent 5f7387c
commit 4d9a3f2
Show file tree

Hide file tree

Showing 3 changed files with 113 additions and 29 deletions.
diff --git a/notebooks/forced_photometry_lightcurve.ipynb b/notebooks/forced_photometry_lightcurve.ipynb
@@ -0,0 +1,83 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "1207192d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import logging\n",
+    "from nuztf.forced_photometry import plot_forced_photometry_lightcurve"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "9fec6764",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "logging.getLogger(\"nuztf\").setLevel(\"DEBUG\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "bcc07f6e",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "INFO:nuztf.cat_match:Using Astropy NED query result for name SDSS J145552.56+414030.4 ([223.96902, 41.67512])\n",
+      "INFO:nuztf.forced_photometry:Requesting IRSA photometry (rerequest is False, cached file does not. This will take some time!\n",
+      "INFO:nuztf.forced_photometry:Lightcurve not found, sleeping for 60.0 seconds. Warning: the lightcurve webpage is only updated hourly.\n",
+      "INFO:nuztf.forced_photometry:Lightcurve not found, sleeping for 60.0 seconds. Warning: the lightcurve webpage is only updated hourly.\n",
+      "INFO:nuztf.forced_photometry:Lightcurve not found, sleeping for 60.0 seconds. Warning: the lightcurve webpage is only updated hourly.\n",
+      "INFO:nuztf.forced_photometry:Lightcurve not found, sleeping for 60.0 seconds. Warning: the lightcurve webpage is only updated hourly.\n"
+     ]
+    }
+   ],
+   "source": [
+    "plot_forced_photometry_lightcurve(\n",
+    "    \"SDSS J145552.56+414030.4\",\n",
+    "    plot_mag=False,\n",
+    "    nu_name=[\"IC220624A\"],\n",
+    "    use_difference_flux=False,\n",
+    "    snr_cut=5.0,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f33ec17f",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "nuztf_env",
+   "language": "python",
+   "name": "nuztf_env"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.12"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/nuztf/forced_photometry.py b/nuztf/forced_photometry.py
@@ -73,7 +73,6 @@ def submit_request(
     req_soup = BeautifulSoup(request.text, "html.parser")
 
     # For one pending job (ie: one request), check forced photometry job tables
-
     req_ra = req_dec = req_jds = req_jde = None
 
     # Grab the table section of the HTML
@@ -96,7 +95,6 @@ def submit_request(
     # Check if request is fulfilled using a request status check and parse with beautifulsoup4
     # Note: Requests older than 30 days will not show up here
     # Iterable, as beautifulsoup can parse html columns and output lists
-
     request_completed = False
 
     output_end_time = output_lc_path = ""
@@ -202,31 +200,45 @@ def load_forced_photometry(
 
     df = pd.read_csv(save_path, comment="#", sep=" ")
 
+    raw_count = len(df)
+
     acceptable_proc = ["0", "62", "63", "255"]
     mask = np.array([str(x) in acceptable_proc for x in df["procstatus"]])
 
+    df = df[mask]
+
+    df = df[df["forcediffimflux,"] != "null"]
+    df = df[
+        (df["infobitssci,"] == 0) & (df["scisigpix,"] < 25) & (df["sciinpseeing,"] < 4)
+    ]
+    df = df[df["dnearestrefsrc,"] < 1]
+
     logger.info(
-        f"Found {np.sum(mask)} entries with procstatus in {acceptable_proc}. "
-        f"Dropping {np.sum(~mask)} other entries."
+        f"Found {len(df)} entries passing quality cuts. "
+        f"Dropping {raw_count-len(df)} other entries."
     )
 
-    df = df[mask]
     df["mjd"] = df["jd,"] - 2400000.5
     df["filtercode"] = [f"z{x[-1].lower()}" for x in df["filter,"]]
 
     if not use_difference_flux:
         # Frank's formulae
-        nearestrefflux = 10.0 ** (0.4 * (df["zpdiff,"] - df["nearestrefmag,"]))
-        nearestreffluxunc = df["nearestrefmagunc,"] * nearestrefflux / 1.0857
+        nearest_ref_flux = 10.0 ** (0.4 * (df["zpdiff,"] - df["nearestrefmag,"]))
 
-        Fluxtot = df["forcediffimflux,"] + nearestrefflux
-        Fluxunctot = np.sqrt(df["forcediffimflux,"] ** 2.0 - nearestreffluxunc**2.0)
-        SNRtot = Fluxtot / Fluxunctot
+        nearest_ref_flux_unc = df["nearestrefmagunc,"] * nearest_ref_flux / 1.0857
 
-        df["mag"] = df["zpdiff,"] - 2.5 * np.log10(Fluxtot)
-        df["magerr"] = 1.0857 / SNRtot
+        flux_tot = df["forcediffimflux,"] + nearest_ref_flux
 
-        snr_mask = SNRtot > snr_cut
+        flux_unc_tot = np.sqrt(
+            df["forcediffimfluxunc,"] ** 2.0 + nearest_ref_flux_unc**2.0
+        )
+
+        snr_tot = flux_tot / flux_unc_tot
+
+        df["mag"] = df["zpdiff,"] - 2.5 * np.log10(flux_tot)
+        df["magerr"] = 1.0857 / snr_tot
+
+        snr_mask = snr_tot > snr_cut
         df = df[snr_mask]
 
         logger.info(df)
@@ -328,14 +340,3 @@ def plot_composite_lightcurve(source_name: str, snr_cut: float = 5.0, **kwargs):
 
     data = Table.from_pandas(df)
     lightcurve_from_science_image(source_name=source_name, data=data, **kwargs)
-
-
-if __name__ == "__main__":
-    logging.getLogger("nuztf").setLevel("DEBUG")
-    plot_forced_photometry_lightcurve(
-        "1RXS J145552.7+414026",
-        plot_mag=True,
-        nu_name=["IC220624A"],
-        use_difference_flux=True,
-        snr_cut=5.0,
-    )
diff --git a/nuztf/plot.py b/nuztf/plot.py
@@ -359,13 +359,13 @@ def lightcurve_from_science_image(
 
     # Start Figure
 
-    plt.figure(figsize=(base_width * 1.1, base_height), dpi=dpi)
+    plt.figure(figsize=(base_width * 1.3, base_height), dpi=dpi)
 
     if expanded_labels:
 
-        ax2 = plt.subplot(111)
+        ax = plt.subplot(111)
 
-        ax = ax2.twiny()
+        ax2 = ax.twiny()
 
     else:
         ax = plt.subplot(111)
@@ -598,7 +598,7 @@ def lightcurve_from_science_image(
         ax2.set_xlim(lmjd, umjd)
 
         if plot_title is not None:
-            ax.set_title(f'ZTF Lightcurve of {plot_title.replace("J", " J")}', y=1.4)
+            ax.set_title(f'ZTF Lightcurve of {plot_title.replace("J", " J")}', y=1.6)
 
         ax2.tick_params(axis="both", which="major", labelsize=big_fontsize)
 
@@ -609,7 +609,7 @@ def lightcurve_from_science_image(
 
     ax.legend(
         loc="upper center",
-        bbox_to_anchor=(0.5, 1.22 + 0.2 * float(expanded_labels)),
+        bbox_to_anchor=(0.5, 1.22 + 0.4 * float(expanded_labels)),
         ncol=3 + len(nu_name),
         fancybox=True,
         fontsize=big_fontsize,