Posterior Skymap Development

fast_response/AlertFollowup.py

@@ -370,4 +370,6 @@ class CascadeFollowup(AlertFollowup):
     _dataset = "GFUOnline_v001p02"
     _fix_index = True
     _float_index = not _fix_index
-    _index = 2.5
+    _index = 2.5
+    _nside = 512
+

fast_response/GWFollowup.py

@@ -351,8 +351,8 @@ def get_best_fit_contour(self, proportions=[0.5,0.9]):
         plt.savefig(self.analysispath + '/' + self.analysisid + 'ts_contours.pdf',bbox_inches='tight', dpi=300)
         plt.close()
 
-    def plot_ontime(self, with_contour=True, contour_files=None, label_events=True):
-        return super().plot_ontime(with_contour=True, contour_files=contour_files, label_events=label_events)
+    def plot_ontime(self, with_contour=True, contour_files=None, label_events=True,custom_events=None):
+        return super().plot_ontime(with_contour=True, contour_files=contour_files, label_events=label_events,custom_events=custom_events)
 
     def write_circular(self):
         """
@@ -541,7 +541,7 @@ def inject_scan(self, ra, dec, ns, poisson=True):
         results = dict([('ts',ts),('ns',ns),('gamma',gamma),('ra',ra),('dec',dec)])
         return (results, events)
 
-    def find_coincident_events(self):
+    def find_coincident_events(self,custom_events=None):
         r"""
         Find "coincident events" for a skymap
         based analysis. These are ALL ontime events,
@@ -573,6 +573,27 @@ def find_coincident_events(self):
             events['ts'][i] = self.ts_scan['TS_spatial_prior_0'][idx[0]]
             events['ns'][i] = self.ts_scan['nsignal'][idx[0]]
             events['gamma'][i] = self.ts_scan['gamma'][idx[0]]
+
+        if(custom_events is not None):
+            '''
+            custom_events.dtype.names = 'time', 'run', 'event', 'ra', 'dec', 'azimuth', 'zenith', 'sigma', 'logE', 'subevent'
+            dtype = [('time',np.float32),('run',np.uint32),('event',np.uint32),('ra',np.float32),
+                ('dec',np.float32),('azimuth',np.float32),('zenith',np.float32),('sigma',np.float32),
+                ('logE',np.float32),('subevent',np.float64)]
+            custom_events=custom_events[['time', 'run', 'event', 'ra', 'dec', 'azimuth', 'zenith', 'sigma', 'logE', 'subevent']]
+            '''
+            custom_events = append_fields(
+                custom_events, names=['in_contour', 'ts', 'ns', 'gamma', 'B'],
+                data=np.empty((5, custom_events['ra'].size)),
+                usemask=False)
+            for i in range(custom_events['ra'].size):
+                custom_events['in_contour'][i]=True
+                custom_events['B'][i] = -1
+                custom_events['ts'][i] = -1
+                custom_events['ns'][i] = -1
+                custom_events['gamma'][i] = -1
+            events=np.concatenate([events,custom_events])
+
 
         self.events_rec_array = events
         self.coincident_events = [dict(zip(events.dtype.names, x)) for x  in events]

fast_response/latex/posterior_report_general.tex

@@ -0,0 +1,74 @@
+\documentclass[titlepage]{article}
+\usepackage{tikz}
+\usepackage{xcolor}
+\usepackage{graphicx}
+\usepackage{grffile}
+\usepackage{longtable}
+\usepackage[margin=1.0in]{geometry}
+
+\setlength{\parindent}{0cm} 
+\include{r_posterior}
+
+\newcommand{\degree}{$^{\circ}$}
+\begin{document}
+
+\begin{titlepage}
+  \centering
+  \vspace{4cm}
+  {\huge\bfseries IceCube Fast-Response \\ Posterior Supplemental Report\par}
+  \vspace{1cm}
+  {\LARGE For Internal Use Only\par}
+  \vfill
+  {\Large Source Name: \\ \itshape\sourcename\par}
+  \vspace{0.5cm}
+  {\Large Observation Date(s):\\ \obsdate \par}
+  \vfill
+  \vspace{1cm} 
+  {\Large Report Generated On:\\ \reportdate \par}
+\end{titlepage}
+
+
+\pagebreak
+\section{Results}
+
+{
+  \centering
+  {\Large Prior Zoomed}
+
+  \includegraphics[width=0.9\textwidth]{\PriorSkymap}
+
+
+}
+\section{Posterior Maps}
+
+{
+  \centering
+  {\Large Posterior Skymap}
+
+  \includegraphics[width=0.9\textwidth]{\PosteriorSkymap}
+  {\Large Zoomed Posterior Skymap}
+
+  \includegraphics[width=0.9\textwidth]{\PosteriorSkymapZoom}
+  {\Large Posterior Skymap W/ Nu Overlay}
+
+  \includegraphics[width=0.9\textwidth]{\PosteriorSkymapZoomWNeutrinos}
+
+
+}
+\pagebreak
+
+{
+  \centering
+  {\Large $N_s$ Profile Plot}
+
+  \includegraphics[width=0.9\textwidth]{\NsProfilePlot}
+
+
+}
+
+
+
+
+\vfill
+
+\end{document}

fast_response/plotting_utils.py

@@ -29,19 +29,19 @@ def plot_zoom(scan, ra, dec, title, reso=3, var="pVal", range=[0, 6],cmap=None):
         cmap = mpl.colors.ListedColormap(pdf_palette)
     hp.gnomview(scan, rot=(np.degrees(ra), np.degrees(dec), 0),
                     cmap=cmap,
+                    cbar=False,
                     max=max(scan),
                     reso=reso,
                     title=title,
                     notext=True,
-                    cbar=False
                     #unit=r""
                     )
 
     plt.plot(4.95/3.*reso*np.radians([-1, 1, 1, -1, -1]), 4.95/3.*reso*np.radians([1, 1, -1, -1, 1]), color="k", ls="-", lw=3)
     hp.graticule(verbose=False)
     plot_labels(dec, ra, reso)
 
-def plot_color_bar(labels=[0.,2.,4.,6.], col_label=r"IceCube Event Time", range=[0,6], cmap=None, offset=-35):
+def plot_color_bar(labels=[0.,2.,4.,6.], col_label=r"IceCube Event Time", range=[0,6], cmap=None, offset=-35,loc=[0.95, 0.2, 0.03, 0.6]):
     """
     Adds a color bar to an existing healpy map
 
@@ -55,9 +55,11 @@ def plot_color_bar(labels=[0.,2.,4.,6.], col_label=r"IceCube Event Time", range=
         colormap to use. if not set default is "Blues"
     offset: int
         offset value for colorbar's label. default is -35
+    loc: float list of length 4
+        location for the colorbar axis in form [x,y,width,height]
     """
     fig = plt.gcf()
-    ax = fig.add_axes([0.95, 0.2, 0.03, 0.6])
+    ax = fig.add_axes(loc)
     labels = labels
     cb = mpl.colorbar.ColorbarBase(ax, cmap="Blues" if cmap is None else cmap,
                         #norm=mpl.colors.Normalize(vmin=range[0], vmax=range[1]), 
@@ -68,7 +70,7 @@ def plot_color_bar(labels=[0.,2.,4.,6.], col_label=r"IceCube Event Time", range=
     cb.set_ticklabels(labels)
     cb.update_ticks()
 
-def plot_labels(src_dec, src_ra, reso):
+def plot_labels(src_dec, src_ra, reso,nPix=200,label_scale=1):
     """
     Add labels to healpy zoom
 
@@ -82,6 +84,7 @@ def plot_labels(src_dec, src_ra, reso):
         Resolution (arcmins)
     """
     fontsize = 20
+    reso=reso*nPix/200
     plt.text(-1*np.radians(1.75*reso),np.radians(0), r"%.1f$^{\circ}$"%(np.degrees(src_dec)),
              horizontalalignment='right',
              verticalalignment='center', fontsize=fontsize)
@@ -100,14 +103,44 @@ def plot_labels(src_dec, src_ra, reso):
     plt.text(np.radians(-reso),np.radians(-1.75*reso), r"%.1f$^{\circ}$"%(reso+np.degrees(src_ra)),
              horizontalalignment='center',
              verticalalignment='top', fontsize=fontsize)
-    plt.text(-1*np.radians(2.35*reso), np.radians(0), r"declination", 
+    plt.text(-1*np.radians(2.35*reso*label_scale), np.radians(0), r"declination", 
                 ha='center', va='center', rotation=90, fontsize=fontsize)
-    plt.text(np.radians(0), np.radians(-2.05*reso), r"right ascension", 
+    plt.text(np.radians(0), np.radians(-2.05*reso*label_scale), r"right ascension", 
                 ha='center', va='center', fontsize=fontsize)
+def compute_ang_err(ra,dec,sigma):
+        dec = np.pi/2 - dec
+        sigma = np.rad2deg(sigma)
+        delta, step, bins = 0, 0, 0
+        delta= sigma/180.0*np.pi
+        step = 1./np.sin(delta)/20.
+        bins = int(360./step)
+        Theta = np.zeros(bins+1, dtype=np.double)
+        Phi = np.zeros(bins+1, dtype=np.double)
+        # define the contour
+        for j in range(0,bins):
+                phi = j*step/180.*np.pi
+                vx = np.cos(phi)*np.sin(ra)*np.sin(delta) + np.cos(ra)*(np.cos(delta)*np.sin(dec) + np.cos(dec)*np.sin(delta)*np.sin(phi))
+                vy = np.cos(delta)*np.sin(dec)*np.sin(ra) + np.sin(delta)*(-np.cos(ra)*np.cos(phi) + np.cos(dec)*np.sin(ra)*np.sin(phi))
+                vz = np.cos(dec)*np.cos(delta) - np.sin(dec)*np.sin(delta)*np.sin(phi)
+                Theta[j], Phi[j] = hp.vec2ang(np.array([vx, vy, vz]))
+
+        Theta[bins] = Theta[0]
+        Phi[bins] = Phi[0]
+
+        return Theta, Phi
+def plot_events2(dec,ra,sigmas, src_ra, src_dec, reso, sigma_scale=5., col = 'k', constant_sigma=False,
+                    same_marker=False, energy_size=False, with_mark=True, with_dash=False,
+                    label='',resolution=0):
+    #plot based on explicitly determined contours rather than using the s parameter. 
+    for i in range(len(ra)):
+        ev_contour = compute_ang_err(ra[i],dec[i],sigmas[i])
+        hp.projplot(ev_contour[0], ev_contour[1], linewidth=1.75, color=col[0],
+                    linestyle="solid",coord='C')
+
 
 def plot_events(dec, ra, sigmas, src_ra, src_dec, reso, sigma_scale=5., col = 'k', constant_sigma=False,
                     same_marker=False, energy_size=False, with_mark=True, with_dash=False,
-                    label=''):
+                    label='',resolution=0.025):
     """
     Adds events to a healpy zoom plot. Events are expected to be from self.llh.exp
 
@@ -146,7 +179,7 @@ def plot_events(dec, ra, sigmas, src_ra, src_dec, reso, sigma_scale=5., col = 'k
 
     if sigma_scale is not None:
         sigma = np.degrees(sigmas)/sigma_scale
-        sizes = 5200*sigma**2
+        sizes = 5200*np.power(resolution/.025,.5)*sigma**2
         if constant_sigma:
             sizes = 20*np.ones_like(sizes)
         if with_dash:
@@ -308,8 +341,8 @@ def make_public_zoom_skymap(skymap, events, ra, dec, with_contour=True, name='te
 
     #plt.scatter(0,0, marker='*', c = 'k', s = 130, label = label_str) 
     #plt.legend(loc = 2, ncol=1, mode = 'expand', fontsize = 18.5, framealpha = 0.95)
+
     plot_color_bar(range=[0,6], cmap=cmap, col_label='GW Map Probability', offset = -10,
-                   labels=[f'{min(skymap):.1e}',f'{max(skymap):.1e}'])
-
+                labels=[f'{min(skymap):.1e}',f'{max(skymap):.1e}'])
     plt.savefig(f'./{name}_skymap_zoom_public.png', bbox_inches='tight', dpi=300)
     plt.close()