Merge 89a5111 into 67519c1

nuztf/base_scanner.py

@@ -18,7 +18,6 @@
 from astropy.coordinates import SkyCoord, Distance
 from astropy.cosmology import FlatLambdaCDM
 from ztfquery import fields as ztfquery_fields
-from ztfquery.fields import FIELD_DATAFRAME
 from gwemopt.ztf_tiling import get_quadrant_ipix
 from ampel.ztf.t0.DecentFilter import DecentFilter
 from ampel.ztf.dev.DevAlertConsumer import DevAlertConsumer
@@ -38,6 +37,7 @@
 from nuztf.plot import lightcurve_from_alert
 from nuztf.utils import cosmo
 from nuztf.fritz import save_source_to_group
+from nuztf.flatpix import get_flatpix
 
 DEBUG = False
 RATELIMIT_CALLS = 10
@@ -61,6 +61,16 @@ def __init__(
         self.cone_nside = cone_nside
         self.t_min = t_min
 
+        (
+            self.map_coords,
+            self.pixel_nos,
+            self.nside,
+            self.map_probs,
+            self.data,
+            self.pixel_area,
+            self.key,
+        ) = self.unpack_skymap()
+
         if not hasattr(self, "prob_threshold"):
             self.prob_threshold = None
 
@@ -121,6 +131,9 @@ def get_name(self):
     def get_full_name(self):
         raise NotImplementedError
 
+    def unpack_skymap(self):
+        raise NotImplementedError
+
     @staticmethod
     def get_tiling_line():
         return ""
@@ -134,7 +147,12 @@ def remove_variability_line():
         raise NotImplementedError
 
     def get_overlap_line(self):
-        raise NotImplementedError
+        """ """
+        return (
+            f"We covered {self.overlap_prob:.1f}% ({self.double_extragalactic_area:.1f} sq deg) "
+            f"of the reported localization region. "
+            "This estimate accounts for chip gaps. "
+        )
 
     def filter_ampel(self, res):
         self.logger.debug("Running AMPEL filter")
@@ -415,7 +433,7 @@ def draft_gcn(self):
             "GROWTH acknowledges generous support of the NSF under PIRE Grant No 1545949.\n"
             "Alert distribution service provided by DIRAC@UW (Patterson et al. 2019).\n"
             "Alert database searches are done by AMPEL (Nordin et al. 2019).\n"
-            "Alert filtering is performed with the AMPEL Follow-up Pipeline (Stein et al. 2021).\n"
+            "Alert filtering is performed with the nuztf (Stein et al. 2021, https://github.com/desy-multimessenger/nuztf).\n"
         )
         if self.dist:
             text += "Alert filtering and follow-up coordination is being undertaken by the Fritz marshal system (FIXME CITATION NEEDED)."
@@ -743,18 +761,7 @@ def __init__(self, data):
         pix_map = dict()
         pix_obs_times = dict()
 
-        infile = os.path.join(
-            "nuztf", "data", f"ztf_fields_ipix_nside={self.nside}.pickle"
-        )
-
-        # Generate a lookup table for field healpix
-        # if none exists (because this is computationally costly)
-        if not os.path.isfile(infile):
-            self.generate_flatpix_file()
-
-        with open(infile, "rb") as f:
-            field_pix = pickle.load(f)
-        f.close()
+        field_pix = get_flatpix(nside=self.nside, logger=self.logger)
 
         for i, obs_time in enumerate(tqdm(obs_times)):
 
@@ -1022,46 +1029,6 @@ def plot_overlap_with_observations(
         )
         return fig, message
 
-    def generate_flatpix_file(self):
-        """
-        Generate and save the fields-healpix lookup table
-        """
-
-        self.logger.info(
-            f"Generating field-healpix lookup table for nside={self.nside}"
-        )
-
-        field_dataframe = FIELD_DATAFRAME.reset_index()
-
-        fields = field_dataframe["ID"].values
-        ras = field_dataframe["RA"].values
-        decs = field_dataframe["Dec"].values
-
-        flat_pix_dict = dict()
-
-        for i, field in tqdm(enumerate(fields), total=len(fields)):
-
-            ra = ras[i]
-            dec = decs[i]
-            pix = get_quadrant_ipix(self.nside, ra, dec)
-
-            flat_pix = []
-
-            for sub_list in pix:
-                for p in sub_list:
-                    flat_pix.append(p)
-
-            flat_pix = list(set(flat_pix))
-            flat_pix_dict[field] = flat_pix
-
-        outdir = os.path.join("nuztf", "data")
-        if not os.path.exists(outdir):
-            os.makedirs(outdir)
-
-        outfile = os.path.join(outdir, f"ztf_fields_ipix_nside={self.nside}.pickle")
-        with open(outfile, "wb") as f:
-            pickle.dump(flat_pix_dict, f)
-
     def crosscheck_prob(self):
 
         try:

nuztf/flatpix.py

@@ -0,0 +1,68 @@
+import logging
+import os
+import pickle
+from tqdm import tqdm
+from ztfquery.fields import FIELD_DATAFRAME
+from gwemopt.ztf_tiling import get_quadrant_ipix
+
+
+def get_flatpix_path(nside: int):
+    outdir = os.path.join("nuztf", "data")
+    if not os.path.exists(outdir):
+        os.makedirs(outdir)
+
+    outfile = os.path.join(outdir, f"ztf_fields_ipix_nside={nside}.pickle")
+    return outfile
+
+
+def generate_flatpix_file(nside: int, logger=logging.getLogger(__name__)):
+    """
+    Generate and save the fields-healpix lookup table
+    """
+
+    logger.info(f"Generating field-healpix lookup table for nside={nside}")
+
+    field_dataframe = FIELD_DATAFRAME.reset_index()
+
+    fields = field_dataframe["ID"].values
+    ras = field_dataframe["RA"].values
+    decs = field_dataframe["Dec"].values
+
+    flat_pix_dict = dict()
+
+    for i, field in tqdm(enumerate(fields), total=len(fields)):
+
+        ra = ras[i]
+        dec = decs[i]
+        pix = get_quadrant_ipix(nside, ra, dec)
+
+        flat_pix = []
+
+        for sub_list in pix:
+            for p in sub_list:
+                flat_pix.append(p)
+
+        flat_pix = list(set(flat_pix))
+        flat_pix_dict[field] = flat_pix
+
+    outfile = get_flatpix_path(nside=nside)
+
+    logger.info(f"Saving to {outfile}")
+
+    with open(outfile, "wb") as f:
+        pickle.dump(flat_pix_dict, f)
+
+
+def get_flatpix(nside: int, logger=logging.getLogger(__name__)):
+
+    infile = get_flatpix_path(nside=nside)
+
+    # Generate a lookup table for field healpix
+    # if none exists (because this is computationally costly)
+    if not os.path.isfile(infile):
+        generate_flatpix_file(nside=nside, logger=logger)
+
+    with open(infile, "rb") as f:
+        field_pix = pickle.load(f)
+
+    return field_pix

nuztf/neutrino_scanner.py

@@ -123,14 +123,6 @@ def __init__(
             * abs(np.cos(np.radians(dec[0])))
         )
         self.logger.info(f"Projected Area: {self.rectangular_area:.3f} sq. deg.")
-        (
-            self.map_coords,
-            self.pixel_nos,
-            self.nside,
-            self.map_probs,
-            self.data,
-            self.key,
-        ) = self.unpack_map()
 
     def get_name(self):
         """ """
@@ -140,13 +132,6 @@ def get_full_name(self):
         """ """
         return f"neutrino event {self.get_name()} ({self.author} et. al, GCN {self.gcn_no})"
 
-    def get_overlap_line(self):
-        """ """
-        return (
-            f"We covered {self.overlap_prob:.1f}% ({self.healpix_area:.1f} sq deg) of the reported localization region. "
-            "This estimate accounts for chip gaps. "
-        )
-
     def candidate_text(
         self, ztf_id: str, first_detection: float, lul_lim: float, lul_jd: float
     ):
@@ -253,9 +238,8 @@ def in_contour(self, ra_deg, dec_deg):
 
         return np.logical_and(in_ra, in_dec)
 
-    def unpack_map(self):
+    def unpack_skymap(self):
         """ """
-        # nside = self.cone_nside
         nside = 1024
         map_coords = []
         pixel_nos = []
@@ -290,4 +274,6 @@ def unpack_map(self):
         data = np.zeros(hp.nside2npix(nside), dtype=np.dtype([(key, float)]))
         data[np.array(pixel_nos)] = map_probs
 
-        return map_coords, pixel_nos, nside, map_probs, data, key
+        pixel_area = hp.nside2pixarea(nside, degrees=True) * float(len(map_coords))
+
+        return map_coords, pixel_nos, nside, map_probs, data, pixel_area, key

nuztf/plot.py

@@ -269,9 +269,6 @@ def absmag_to_mag(absmag):
 
     # Ugly hack because secondary_axis does not work with astropy.time.Time datetime conversion
     mjd_min = min(np.min(df.mjd.values), t_0_mjd)
-
-    print(mjd_min, t_0_mjd)
-
     mjd_max = max(np.max(df.mjd.values), t_0_mjd)
     length = mjd_max - mjd_min
 

nuztf/skymap_scanner.py

@@ -142,13 +142,6 @@ def __init__(
         self.logger.info("Reading map")
 
         self.pixel_threshold = self.find_pixel_threshold(self.data[self.key])
-        (
-            self.map_coords,
-            self.pixel_nos,
-            self.map_probs,
-            self.nside,
-            self.pixel_area,
-        ) = self.unpack_skymap()
 
         BaseScanner.__init__(
             self,
@@ -168,6 +161,12 @@ def __init__(
         if not os.path.exists(self.cache_dir):
             os.makedirs(self.cache_dir)
 
+    def get_full_name(self):
+        if self.event_name is not None:
+            return self.event_name
+        else:
+            return "?????"
+
     def get_alerts(self):
         """Scan the skymap area and get ZTF transients"""
         self.logger.info("Commencing skymap scan")
@@ -320,16 +319,6 @@ def get_obs_line(self):
         """ """
         return "Each exposure was 30s with a typical depth of 20.5 mag."
 
-    def get_overlap_line(self):
-        """ """
-        return (
-            "We covered {0:.1f}% of the enclosed probability "
-            "based on the map in {1:.1f} sq deg. "
-            "This estimate accounts for chip gaps. ".format(
-                self.overlap_prob, self.area
-            )
-        )
-
     @staticmethod
     def remove_variability_line():
         """ """
@@ -702,7 +691,15 @@ def unpack_skymap(self):
             map_coords, dtype=np.dtype([("ra", float), ("dec", float)])
         )
 
-        return map_coords, pixel_nos, self.data[self.key][mask], nside, pixel_area
+        return (
+            map_coords,
+            pixel_nos,
+            nside,
+            self.data[self.key][mask],
+            self.data,
+            pixel_area,
+            self.key,
+        )
 
     def find_cone_coords(self):
         """ """

tests/test_neutrino_scanner.py

@@ -113,9 +113,7 @@ def test_scan(self):
         print(repr(res))
 
         # Update the true using repr(res)
-        true_gcn = f"Astronomer Name (Institute of Somewhere), ............. report,\n\nOn behalf of the Zwicky Transient Facility (ZTF) and Global Relay of Observatories Watching Transients Happen (GROWTH) collaborations: \n\nAs part of the ZTF neutrino follow up program (Stein et al. 2022), we observed the localization region of the neutrino event IceCube-200620A (Santander et. al, GCN 27997) with the Palomar 48-inch telescope, equipped with the 47 square degree ZTF camera (Bellm et al. 2019, Graham et al. 2019). We started observations in the g- and r-band beginning at 2020-06-21 04:53 UTC, approximately 25.8 hours after event time. We covered 77.7% (1.2 sq deg) of the reported localization region. This estimate accounts for chip gaps. Each exposure was 300s with a typical depth of 21.0 mag. \n \nThe images were processed in real-time through the ZTF reduction and image subtraction pipelines at IPAC to search for potential counterparts (Masci et al. 2019). AMPEL (Nordin et al. 2019, Stein et al. 2021) was used to search the alerts database for candidates. We reject stellar sources (Tachibana and Miller 2018) and moving objects, and apply machine learning algorithms (Mahabal et al. 2019) . We are left with the following high-significance transient candidates by our pipeline, all lying within the 90.0% localization of the skymap.\n\n+--------------------------------------------------------------------------------+\n| ZTF Name     | IAU Name  | RA (deg)    | DEC (deg)   | Filter | Mag   | MagErr |\n+--------------------------------------------------------------------------------+\n| ZTF18acvhwtf | AT2020ncs | {hist_and_new_values['ZTF18acvhwtf']['ra']:011.7f} | {hist_and_new_values['ZTF18acvhwtf']['dec']:+011.7f} | r      | {hist_and_new_values['ZTF18acvhwtf']['mag']:.2f} | {hist_and_new_values['ZTF18acvhwtf']['mag_err']:.2f}   | {old_flag} \n| ZTF20abgvabi | AT2020ncr | 162.5306341 | +12.1461187 | g      | 20.58 | 0.19   | (MORE THAN ONE DAY SINCE SECOND DETECTION) \n+--------------------------------------------------------------------------------+\n\n \n\nAmongst our candidates, \n\nZTF18acvhwtf was first detected on 2458461.9815278. It has a spec-z of 0.291 [{hist_and_new_values['ZTF18acvhwtf']['z_dist']:.0f} Mpc] and an abs. mag of {hist_and_new_values['ZTF18acvhwtf']['absmag']:.1f}. Distance to SDSS galaxy is {hist_and_new_values['ZTF18acvhwtf']['ned_dist_new']:.2f} arcsec. [MILLIQUAS: SDSS J104816.25+120734.7 - 'Q'-type source ({hist_and_new_values['ZTF18acvhwtf']['milliquas_dist_new']:.2f} arsec)] [TNS NAME=AT2020ncs]\nZTF20abgvabi was first detected on 2458995.6705903. WISEA J105007.28+120846.1 ['G'-type source (0.00 arsec)] [TNS NAME=AT2020ncr]\n\n\nZTF and GROWTH are worldwide collaborations comprising Caltech, USA; IPAC, USA; WIS, Israel; OKC, Sweden; JSI/UMd, USA; DESY, Germany; TANGO, Taiwan; UW Milwaukee, USA; LANL, USA; TCD, Ireland; IN2P3, France.\n\nGROWTH acknowledges generous support of the NSF under PIRE Grant No 1545949.\nAlert distribution service provided by DIRAC@UW (Patterson et al. 2019).\nAlert database searches are done by AMPEL (Nordin et al. 2019).\nAlert filtering is performed with the AMPEL Follow-up Pipeline (Stein et al. 2021).\n"
-
-        print(res)
+        true_gcn = f"Astronomer Name (Institute of Somewhere), ............. report,\n\nOn behalf of the Zwicky Transient Facility (ZTF) and Global Relay of Observatories Watching Transients Happen (GROWTH) collaborations: \n\nAs part of the ZTF neutrino follow up program (Stein et al. 2022), we observed the localization region of the neutrino event IceCube-200620A (Santander et. al, GCN 27997) with the Palomar 48-inch telescope, equipped with the 47 square degree ZTF camera (Bellm et al. 2019, Graham et al. 2019). We started observations in the g- and r-band beginning at 2020-06-21 04:53 UTC, approximately 25.8 hours after event time. We covered 77.7% (1.2 sq deg) of the reported localization region. This estimate accounts for chip gaps. Each exposure was 300s with a typical depth of 21.0 mag. \n \nThe images were processed in real-time through the ZTF reduction and image subtraction pipelines at IPAC to search for potential counterparts (Masci et al. 2019). AMPEL (Nordin et al. 2019, Stein et al. 2021) was used to search the alerts database for candidates. We reject stellar sources (Tachibana and Miller 2018) and moving objects, and apply machine learning algorithms (Mahabal et al. 2019) . We are left with the following high-significance transient candidates by our pipeline, all lying within the 90.0% localization of the skymap.\n\n+--------------------------------------------------------------------------------+\n| ZTF Name     | IAU Name  | RA (deg)    | DEC (deg)   | Filter | Mag   | MagErr |\n+--------------------------------------------------------------------------------+\n| ZTF18acvhwtf | AT2020ncs | {hist_and_new_values['ZTF18acvhwtf']['ra']:011.7f} | {hist_and_new_values['ZTF18acvhwtf']['dec']:+011.7f} | r      | {hist_and_new_values['ZTF18acvhwtf']['mag']:.2f} | {hist_and_new_values['ZTF18acvhwtf']['mag_err']:.2f}   | {old_flag} \n| ZTF20abgvabi | AT2020ncr | 162.5306341 | +12.1461187 | g      | 20.58 | 0.19   | (MORE THAN ONE DAY SINCE SECOND DETECTION) \n+--------------------------------------------------------------------------------+\n\n \n\nAmongst our candidates, \n\nZTF18acvhwtf was first detected on 2458461.9815278. It has a spec-z of 0.291 [{hist_and_new_values['ZTF18acvhwtf']['z_dist']:.0f} Mpc] and an abs. mag of {hist_and_new_values['ZTF18acvhwtf']['absmag']:.1f}. Distance to SDSS galaxy is {hist_and_new_values['ZTF18acvhwtf']['ned_dist_new']:.2f} arcsec. [MILLIQUAS: SDSS J104816.25+120734.7 - 'Q'-type source ({hist_and_new_values['ZTF18acvhwtf']['milliquas_dist_new']:.2f} arsec)] [TNS NAME=AT2020ncs]\nZTF20abgvabi was first detected on 2458995.6705903. WISEA J105007.28+120846.1 ['G'-type source (0.00 arsec)] [TNS NAME=AT2020ncr]\n\n\nZTF and GROWTH are worldwide collaborations comprising Caltech, USA; IPAC, USA; WIS, Israel; OKC, Sweden; JSI/UMd, USA; DESY, Germany; TANGO, Taiwan; UW Milwaukee, USA; LANL, USA; TCD, Ireland; IN2P3, France.\n\nGROWTH acknowledges generous support of the NSF under PIRE Grant No 1545949.\nAlert distribution service provided by DIRAC@UW (Patterson et al. 2019).\nAlert database searches are done by AMPEL (Nordin et al. 2019).\nAlert filtering is performed with the nuztf (Stein et al. 2021, https://github.com/desy-multimessenger/nuztf).\n"
 
         self.assertEqual(res, true_gcn)